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1 ABSTRACT  (CONTINUE  ON  REVERSE  SIDE  IF  NECESSARY  AND  IDENTIFY  BY  BLOCK  NUMBER) 

An  advanced  version  of  the  EPIC-3  computer  program  is  described.  The  numerical 
technique  is  based  on  a Lagrangian  finite  element  formulation  where  the  equations  of 
motion  are  integrated  directly  rather  than  through  the  traditional  stiffness  matrix 
approach.  It  is  primarily  intended  for  three  dimensional  analysis  of  penetration  and 
wave  propagation  problems  resulting  from  high  velocity  impact,  although  other 
applications  are  possible.  The  two  most  significant  improvements  over  the  initial 
version  of  the  EPIC-3  code,  are  that  problems  of  unlimited  size  can  be  run  and  that  _ 
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20.  ABSTRACT  (Continued) 

the  code  is  programmed  to  eventually  run  on  fourth  generation  vector  computers. 
Other  improvements  include  the  capability  to  apply  external  pressures,  reduction  of 
required  input,  expanded  geometry  generators,  a more  general  sliding  surface 
capability  and  expanded  plotting  capability. 
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This  final  report  on  the  development  of  the  advanced  EPIC-3  code  was 
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for  their  assistance  in  running  the  sample  problems. 

This  report  has  been  reviewed  by  the  Information  Officer  (01)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS) . At  NTIS 
it  will  be  available  to  the  general  public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 
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SECTION  I 

INTRODUCTION  AND  SUMMARY 

This  report  documents  an  advanced  version  of  the  EPIC-3  computer  pro- 
gram (Elastic-Plastic  hnpact  Computations  in  3^  dimensions).  It  is  pri- 
marily intended  for  analysis  of  three-dimensional  wave  propagation  and 
penetration  problems  resulting  from  high-velocity  impact.  It  is  based  on 
an  explicit  Lagrangian  finite-element  formulation  where  the  equations  of 
motion  are  integrated  directly  rather  than  through  the  traditional  stiffness 
matrix  approach.  This  code  is  an  extension  of  an  earlier  version  of  EPIC- 
3 (Reference  1)  with  significant  improvements  which  are  summarized  later 
in  this  section. 

The  EPIC-3  code  has  material  descriptions  which  include  strain  hardening, 
strain  rate  effects,  thermal  softening  and  fracture.  Geometry  generators 
are  included  to  quickly  generate  flat  plates,  spheres  and  rods  with  various 
nose  shapes.  It  has  the  capability  to  include  sliding  surfaces,  and  it  pro- 
vides plots  of  deformed  geometry  and  velocity  vectors.  In  addition,  the 
finite-element  computational  approach  has  several  specific  advantages  when 
compared  to  the  more  commonly  used  finite-difference  approach: 

• The  finite-element  formulation  does  not  require  an  orderly 
grid  and  is  therefore  well-suited  to  represent  complicated 
geometrical  shapes  by  use  of  three-dimensional  tetrahedron 
elements. 

• The  finite-element  formulation  allows  the  use  of  tetrahedron 
elements  which  are  well-suited  to  represent  severe  distor- 
tions. 
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• The  tetrahedron  elements  are  in  a state  of  constant  strain 
such  that  all  material  in  an  element  behaves  uniformly. 

This  allows  for  an  accurate  and  convenient  selection  of  con- 
stant stresses  and  pressures  within  the  element. 

• The  finite-element  formulation  is  well-suited  for  various 
boundary  conditions.  Free  boundaries  require  no  special 
accommodations,  and  restrained  boundaries  are  simply 
represented  by  setting  the  appropriate  nodal  displacements 
and  velocities  equal  to  zero. 

Although  the  advanced  EPIC-3  code  has  many  similarities  to  the  earlier 
version  (Reference  1),  there  are  many  specific  extensions  and  improve- 
ments. A summary  of  these  follows: 

• The  better-structured  code  consists  of  three  separate  pro- 
grams: Preprocessor,  Main  Routine,  and  Postprocessor. 
These  three  programs  are  connected  by  a common  restart 
tape. 

• Provided  the  nodal  bandwidth  can  be  core  contained,  prob- 
lems of  unlimited  size  can  be  run  by  buffering  data  between 
disk  files  and  core.  For  smaller  problems,  where  all  nodes 
can  be  core  contained,  an  option  is  provided  to  bypass  this 
buffering. 

• The  code  is  programmed  to  eventually  run  on  fourth-gener- 
ation vector  computers.  It  is  expected  that  it  will  also  run 
faster  on  the  CDC  6600  and  7600  computers  since  these 
machines  also  have  some  vectorization  characteristics. 

• Capability  is  provided  to  apply  external  pressures  on  the 
triangular  faces  of  the  tetrahedron  elements.  Pressures 
can  also  be  applied  as  a function  of  time. 
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• Elastic  snapback  is  included  to  allow  slower  problems 

to  be  run  more  accurately.  Internal  energy  computations 
are  also  improved  to  better  conserve  total  energy  (kinetic 
plus  internal)  under  all  conditions. 

• Reduced  input  is  required  for  sliding  surfaces  and  geometry 
generators. 

• Expanded  geometry  generators  include  solid  and  hollow  rods 
with  various  nose  shapes,  solid  and  hollow  spheres,  and  flat 
plates  with  variable  expansion  factors.  There  are  no  limita- 
tions on  the  number  of  rings  in  the  rods,  noses  and  spheres. 

• Expanded  sliding-surface  capability  includes  frictional 
effects.  A general  search  routine  is  also  provided  for  cases 
where  the  top  surface  of  the  target  is  not  a single- valued 
function  along  lines  normal  to  the  top  of  the  undeformed 
target. 

• Expanded  plotting  capability  provides  two-  and  three-dimen- 
sional plots  of  geometry  and  velocity  vectors. 

The  following  sections  of  this  report  present  the  formulation,  a description 
of  the  computer  program,  user  instructions,  and  example  problems. 
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SECTION  II 
FORMULATION 


1.  COMPUTATIONAL  TECHNIQUE 


The  EPIC-3  computational  technique  is  shown  schematically  in  Figure  1. 

The  first  step  in  the  process  is  to  represent  the  geometry  with  tetrahedron 
elements  having  specific  material  characteristics.  Then  the  distributed 
mass  is  lumped  at  the  nodes  (element  corners),  and  initial  velocities  are 
assigned  to  represent  the  motion  at  impact.  If  the  problem  involves  applied 
pressures,  there  may  be  no  initial  velocities,  and  the  initial  conditions  are 
established  by  the  application  of  the  pressures. 


REPRESENT  GEOMETRY 
WITH  ELEMENTS 


LUMP  MASS 
AT  NODES 


STRAINS  & 
STRAIN  RATES 


STRESSES 


L J 

ASSIGN 

u 

I—  

DISPLACEMENTS 

— B 

INITIAL  CONDITIONS 


INTEGRATION  LOOP 


Figure  1.  EPIC- 3 Computational  Technique 
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After  the  initial  conditions  are  established,  the  integration  loop  begins  as 
shown  in  Figure  1.  The  first  step  is  to  obtain  displacements  and  velocities 
of  the  nodes.  If  it  is  assumed  the  lines  connecting  the  nodes  (element  edges) 
remain  straight,  then  the  displacements  and  velocities  within  the  elements 
must  vary  linearly.  From  these  displacements  and  velocities,  the  strains 
and  strain  rates  within  the  elements  can  be  obtained.  Since  the  strains 
and  strain  rates  are  derivatives  of  Unear  displacement  and  velocity  functions 
within  the  elements,  the  resulting  strains  and  strain  rates  are  constant  with- 
in the  elements. 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates, 
internal  energies  and  material  properties.  Since  the  strains  and  strain 
rates  are  constant  within  the  elements,  the  stresses  are  also  constant.  The 
stresses  are  obtained  by  combining  elastic  or  plastic  deviator  stresses  with 
hydrostatic  pressure.  The  deviator  stresses  represent  the  shear  strength 
capability  of  the  material,  and  the  hydrostatic  pressure  is  obtained  from  the 
volumetric  strain  and  internal  energy  of  the  element.  An  artificial  viscosity 
is  also  included  to  damp  out  localized  oscillations  caused  by  representing 
continuous  media  with  lumped  masses. 

After  the  element  stresses  are  determined,  it  is  necessary  to  obtain  con- 
centrated forces  at  the  nodes.  These  forces  are  statically  equivalent  to 
the  distributed  stresses  within  the  element  and  the  applied  pressures  acting 
on  the  faces  of  the  element.  They  are  dependent  on  the  element  geometry 
and  the  magnitude  of  the  stresses  and  pressures.  When  the  concentrated 
forces  are  applied  to  the  concentrated  masses,  the  nodal  accelerations  are 
defined,  and  the  equations  of  motion  are  applied  to  determine  new  displace- 
ments and  velocities.  The  integration  loop  is  then  repeated  until  the  time 
of  interest  has  elapsed. 
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Another  feature  of  the  basic  technique  is  the  ability  to  represent  sliding 
between  two  surfaces.  This  is  accomplished  with  a momentum  exchange 
principle  which  allows  for  closing,  sliding,  and  separation  of  the  two  sur- 
faces. It  should  be  noted  that  the  integration  time  increment  must  be 
properly  controlled  to  prevent  numerical  instability.  This  is  accomplished 
by  limiting  the  time  increment  to  a fraction  of  the  time  required  to  travel 
across  the  minimum  altitude  of  the  element  at  the  sound  velocity  of  the 
material. 


2.  ELEMENT  GEOMETRY 

A typical  tetrahedron  element  is  shown  in  Figure  2.  It  is  geometrically 

defined  by  nodes  i,  j,  m,  and  p.  The  formulation  is  based  on  nodes  i,  j, 

and  m being  positioned  in  a counterclockwise  manner  when  viewed  from 

node  p.  The  mass  of  each  element  is  V p where  V and  p are  the  initial 

o o o o 

volume  and  density  of  the  element.  This  mass  is  equally  distributed  to  con- 
centrated masses  at  the  four  nodes  such  that  the  total  mass  at  node  i,  M., 

1 

is  equal  to  one-fourth  the  mass  of  all  elements  which  contain  that  node. 


The  coordinates  of  node  i are  designated  x.,  y^,  z.,  and  the  corresponding 

velocities  are  designated  u.,  v.,  w.. 
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3.  STRAINS  AND  STRAIN  RATES 


The  incremental  strains  which  occur  during  each  cycle  of  integration  are 
obtained  by  multiplying  the  strain  rates  by  the  integration  time  increment. 
The  strain  rates  are  obtained  from  the  current  geometry  of  the  element 
and  the  velocities  of  the  nodes.  If  it  is  assumed  the  velocities  vary  linearly 
between  the  nodes,  the  x,  y,  and  z velocities  (u,  v,  w)  within  each  element 
can  be  expressed  as 


u = + (*2X  + <*gY  + a4z 


v = + c*gX  + a^y  + OgZ 


w = «q  + «10x  +Q,iiY  + 


11J  12 


(1) 

(2) 

(3) 


where  . . . a ^ are  geometry-  and  velocity-dependent  constants  for  each 
element.  It  is  possible  to  solve  for  ...  by  substituting  the  x velocities 
and  coordinates  of  the  four  nodes  into  Equation  (1).  This  gives  four  equations 
and  four  unknowns  so  that  the  constants  (a^  . . . a can  be  evaluated. 

Equation  (1)  can  then  be  expressed  in  terms  of  element  geometry  and  nodal 
velocities. 


u = 


— [(aA  +b.x  + Ciy  + diZ)  u. 


6V 

+ 

+ 

+ 


(a.  + b.x  + c.y  + d.z)  u . 

3 3 3 3 3 

(a  +b  x + c y + d z)u 
m m m m 

(a  + bx  + cy  + dz)u  1 
P P P^  p ' p I 


(4) 
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where  the  volume  is 


and  the  other  geometry-dependent  constants  are 


(6a) 


(6b) 


(6c) 


(6d) 
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The  remaining  geometry-dependent  constants  (a.,  b.,  c.,  d.,  etc.  ) are  ob- 

3 J J J 

tained  by  a systematic  interchange  of  signs  and  subscripts.  The  y and  z 
velocities  in  Equations  (2)  and  (3)  are  obtained  in  a similar  manner  and  are 
identical  to  Equation  (1)  except  the  x velocities  at  the  four  nodes  are  re- 
placed by  the  y and  z velocities. 


After  the  velocities  are  obtained,  it  is  possible  to  determine  the  normal 

strain  rates  (e  , c , e ),  the  shear  strain  rates  (y  , y , y ) and  the 
x y z xy  xz’  'yz 

spin  rates  (id  . u , u>  ) of  the  element, 
x y z 


i =|a 

(7) 

• d v 

ey  = 

(8) 

• 3w 

C*  ” 1 — ■ 

"z  3z 

(9) 

3u  By 
yxy  ” By  3x 

(10) 

_ 3u  + 3w 
yxz  3z  3x 

(ID 

• • 

• 3v  , 3w 

yz  3z  By 

(12) 

• • 

_ 1 /3w  3v\ 
ux  2 V By  3z  / 

(13) 

1 / 3u  Bw\ 
uy  = 2'3z  " 3x/ 

(14) 

• • 

1 / 3v  3u \ 
wz  2 ' 3x  * By  / 

(15) 

It  can  be  seen  that  Equations  (7)  through  (15)  are  derivatives  of  linear  func- 
tions and  therefore  give  constant  values  within  the  element.  The  volumetric 
strain  is  also  constant  within  the  element  and  is  expressed  as  = V/V  -1, 
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where  V and  Vq  represent  the  current  and  initial  volumes  of  the  element. 

It  should  also  be  noted  that  subsequent  computations  will  involve  deviator 
strain  rates  (e  , e , e ) which  are  readily  obtained  from  the  normal  strain 

. . x.  y z 

rates  (e  , e , e ) by  subtracting  the  average  normal  strain  rate, 
x y z 


4.  STRESSES  AND  PRESSURES 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates, 

internal  energies  and  material  properties  (Reference  2).  The  three  normal 

stresses  (a  , a , a ) are  expressed  in  terms  of  deviator  stresses  (s  , s , 
x y z x y 

sz),  hydrostatic  pressure,  P,  and  artificial  viscosity,  Q: 


°x  = Sx~  (P  +Q) 


ay  = sy-(P+Q) 


az  = Sz  ' (P  +Q) 


(16) 

(17) 

(18) 


Trial  values  of  the  deviator  stresses  and  shear  stresses  at  time  = t + At 
are 


gt  + At  = gt  + 2G(;  A(;  + As 

X XXX 


st  + At  = s1  + 2Ge  At  + As 

y y y y 

sfc  + At  _ at  + 2Ge  At  + As 
z z z z 


t + At  t . „ • . 

T = T + Gy  At  + At 
xy  xy  xy  xy 


(19) 

(20) 
(21) 
(22) 


(23) 


t + 

At  t 

= T 

+ Gy 

At  + At 

xz 

xz 

' xz 

xz 

t + 

At  t 

- T „ 

+ Gy  „ 

At  + At 

yz 

yz 

yz 

yz 

(24) 


In  Equation  (19)  the  first  term  (s^)  is  the  normal  stress  at  the  previous  time 
and  the  second  term  (2GexAt)  is  the  incremental  stress  due  to  the  incremen- 
tal strain  (exAt)  during  that  time  increment,  where  G is  the  elastic  shear 
modulus  and  At  is  the  integration  time  increment.  The  third  term  (Asx)  is 
due  to  shear  stresses  from  the  previous  time  increment,  which  now  act  as 
normal  stresses  due  to  the  new  orientation  of  the  element  caused  by  an  in- 
cremental rotation  (u  At,  iozAt)  during  the  time  increment.  The  remaining 
normal  stresses  and  shear  stresses  have  a similar  form. 


The  correction  terms  for  element  rotations  are 


AS*  ' 2 <“/»  - Viy>  St 

i3y  ’ 2 <“zTxy  - Vyz'1  “ 


As  = 2 (w  t*"  - u Tt  ) At 

z x yz  y xz 

At  =[u  (s*'-s^)+wt*:  - u t1  ] At 

xy  7 x y y yz  x xz 

At  = [u  (s^  - s^)  + u t^  - u t^  ] At 

xz  y z x x xy  z yz 

At  = [oj  (sl  - sl)  +u  Tl  - u Tl  ] At 

yz  x y z z xz  y xy 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 


It  should  be  recalled  that  Equations  (19)  through  (24)  represent  trial  values 
of  the  stresses,  and  they  may  need  to  be  reduced  if  they  violate  the  Von 
Mises  yield  criterion.  An  equivalent  stress  is  given  by 
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3- V 


|*(82  + S2  + S2)  + 3 (T2  +T2  +T2) 

2 x y z xy  xz  yz 


(31) 


If  a is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  S, 
(which  is  discussed  Later  in  this  section),  the  final  deviator  and  shear 
stresses  are  as  given  in  Equations  (19)  through  (24),  If  a is  greater  than 
S,  then  the  stresses  in  Equations  (19)  through  (24)  should  be  multiplied  by 
the  factor  (S/ct).  When  the  reduced  deviator  and  shear  stresses  are  put  in- 
to Equation  (31),  the  result  is  always  a = S, 

During  plastic  flow  it  is  sometimes  necessary  to  determine  the  equivalent 
plastic  strain  for  strain  hardening  effects  on  the  strength  of  the  material 
or  to  determine  if  the  material  has  failed.  The  first  step  in  this  process 
is  to  adjust  the  strain  rates  to  plastic  strain  rates  by  substracting  out  the 
elastic  portion  of  the  strain  rates: 


t + At  t 
s - s - As 

■p  • X X X 

ex  ex  2GAt 


(32) 


(37) 


t + At  t 
t „ - t „ - At 

yz  yz  yz 

crar 


The  equivalent  plastic  strain  rate  is  expressed  as 


C(e  P - e P)  + (ep-ef)  + (ep-~p 


e H) 

y z 


(38) 


+-§<;p  2 + ;p  2 + ;p  2>: 

2 xy  T xz  'yz 


The  equivalent  plastic  strain,  e^,  is  then  obtained  by  integrating  e with  re- 
spect to  time: 


7 + At  = 7 fc  + e At 

P P P 


(39) 


The  equivalent  tensile  strength  of  the  material  may  be  dependent  on  many 
factors,  including  strain,  strain  rate,  pressure  and  temperature.  It  is 
well  known  that  many  materials  behave  differently  under  dynamic  impact 
than  under  static  testing  conditions.  With  few  exceptions,  however,  a pre- 
cise definition  of  material  behavior  under  dynamic  conditions  is  not  avail- 
able. There  is  also  much  to  be  learned  about  fracture  characteristics 
under  these  dynamic  conditions.  The  current  version  of  EPIC-3  allows  the 
equivalent  tensile  strength  to  be  determined  from 


S = S-  [1  + Cj  log  (e)]  [1  + C2P  + C3P2]  [C4  + C5T]  (40) 

p 

In  this  formulation,  S^r  is  generally  taken  to  be  the  statLC  stress,  which  is 
dependent  on  the  equivalent  plastic  strain  of  Equation  (39).  The  three  brack- 
eted terms  allow  the  static  stress  to  be  altered,  based  on  strain  rate, 
pressure  and  temperature.  If  the  constants  are  = C,.  = 0 and 

= 1.  0,  then  S = S^,  which  is  the  strain-dependent  static  stress.  The 
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first  bracketed  term  can  increase  the  strength  due  to  the  equivalent  strain 
rate,  e,  which  is  similar  to  that  described  in  Equation  (38),  but  includes 
the  elastic  portion  of  the  strain  rates.  In  the  second  bracketed  term,  P = 

P + Q,  so  the  effect  of  this  term  is  to  increase  the  strength  due  to  pressure. 
The  final  bracketed  term  includes  the  temperature,  T,  of  the  element, 
which  is  obtained  from 


T = To  + cT  <41) 

s o 

where  Tq  is  the  initial  temperature  of  the  element,  Eg  is  the  internal 

energy  per  original  unit  volume  (defined  later),  C is  the  specific  heat  and 

s 

pQ  is  the  initial  density. 

The  material  is  assumed  to  fracture  if  a specified  value  of  the  equivalent 
plastic  strain,  (Equation  39)  or  the  volumetric  strain,  e ^ - V/VQ-1,  is 
exceeded.  This  is  a very  approximate  failure  model  which  should  be  ex- 
panded when  better  models  become  available.  When  the  failure  criterion 
has  been  met  (exceeded  value  of  or  ^ ),  the  equivalent  tensile  stress  is 
set  to  zero  so  no  shear  stresses  can  be  developed  in  the  failed  element. 
Likewise,  no  tensile  stresses  are  allowed  to  develop.  The  net  result  is 
that  a failed  element  tends  to  act  like  a liquid  inasmuch  as  it  can  develop 
only  hydrostatic  compression  with  no  shear  or  tensile  stresses.  Another 
option  is  also  available  which  sets  all  stresses  (including  the  pressure) 
equal  to  zero. 

The  hydrostatic  pressure  is  dependent  on  the  volumetric  strain  and  the  in- 
ternal energy  in  the  element.  The  EPIC-3  code  uses  the  Mie-Griineisen 
equation  of  state  (Reference  3)  in  the  general  form  P = Pv  + TEg  (1  + p), 
with  the  complete  expression  given  by 


P = (Ktp  + K2p2  + K3h3)  (1  - -t)  + TEg  (1  + p) 


1 4 


where  p = V /V-l,  K^,  K^,  and  are  material-dependent  constants  and 
F is  the  Griineisen  coefficient. 

Since  the  pressure  can  be  significantly  affected  by  the  internal  energy  E , 

s 

it  is  desirable  to  solve  the  pressure  and  energy  equations  simultaneously. 
This  gives 


t + At  Es  - • 5 £<P  + Q)1  + Qfc  + At  + K+  At]  "vM  + AEd 

E = s JE 1 (43) 

s 

1 + . 5T  (1  + p)  e^At 


where  ev  is  the  volumetric  strain  rate  and  AE^  is  the  internal  energy  gen- 
erated by  the  deviator  and  shear  stresses  during  the  previous  cycle. 


AE,  = (i  e +se  +se  +t  y +t  y + r y ) 

d xx  y y z z xy' xy  xz'xz  yzryz 


(V/V0)  At 


(44) 


The  bars  on  the  deviator  and  shear  stresses  and  the  volume  represent 
averages  of  these  values  at  times  t and  t + At.  After  the  total  internal 
energy  at  time  t + At  has  been  determined  from  Equation  (43),  the  pressure 
at  time  t + At  is  determined  from  Equation  (42). 


The  artificial  viscosity  is  combined  with  the  normal  stresses  to  damp  out 
localized  oscillations  of  the  concentrated  masses.  It  tends  to  eliminate 
spurious  oscillations  which  would  otherwise  occur  for  wave  propagation 
problems.  This  technique  was  originally  proposed  by  Von  Neumann  and 
Richtmyer  (Reference  4)  and  has  been  expanded  for  use  in  various  computer 
codes  (Reference  5).  It  is  expressed  in  terms  of  linear  and  quadratic  com- 
ponents and  is  applied  only  when  the  volumetric  strain  rate  is  negative. 
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for  g < 0 
v 


(45) 


W 


Q = Cj  pcgh  livl  +CQph2Uv>2 

Q = 0 for  e ^ 0 
v 

where  cg  is  the  sound  velocity  of  the  material  and  h is  the  minimum  altitude 
of  the  tetrahedron.  Cj  and  Cq  are  the  linear  and  quadratic  dimensionless 
coefficients. 

The  minimum  altitude  of  the  tetrahedron  can  readily  be  obtained  from  the 
previously  calculated  geometry-dependent  constants.  The  altitude  from 
node  i to  the  plane  defined  by  the  other  three  nodes  of  the  tetrahedron  is 


6 V 


, 2 ,2 
+ c . + d . 

l l 


(46) 


The  other  altitudes  (In,  hm>  h^)  are  obtained  by  appropriately  changing  the 
subscripts  of  the  geometry-dependent  constants. 


The  sound  velocities  (Reference  6)  are  obtained  from 


= ^-[Kj  (1  - Tp)  + K2  (2p  - 1.  5fp2) 


(47) 


+ K3  (3p2  - 2 Tp3)  + rEg  + rp/(l  + p)] 


The  sound  velocities  are  also  used  for  determination  of  the  integration  time 
increment. 


1 6 


5.  CONCENTRATED  FORCES 


After  the  element  stresses  are  determined,  the  concentrated  nodal  forces 
can  be  obtained.  These  forces  are  statically  equivalent  to  the  distributed 
stresses  within  the  element  and  the  applied  pressures  acting  on  the  faces  of 
the  element.  They  are  dependent  on  the  displaced  element  geometry  and 
the  magnitude  of  the  stresses  and  pressures.  The  forces  in  the  x,  y and  z 
directions  at  node  i of  an  element  are  given  by 


Fl  =4-[F.  (b.P . + b P + b P ) - (b.o  + c.T  + d.T  j]  (48) 

x 6 t j j mm  ppix  i xy  ixz 

F1  =4-[Ft  (c.P.  + c P +cP)-(c.cr  +b.T  + d.T  )3  (49) 

y 6 t j J mm  p p'  i y i xy  i yz 


Fz  4 CFt  (djpi  + dmPm  + dpPp>  ' (dids  + CiV  + blT»)3  (50) 


The  geometry-dependent  constants  (b.,  c.  d^)  are  again  identical  to  those 
used  for  calculation  of  the  strain  rates  and  altitudes.  The  terms  in  the 


first  sets  of  parentheses  represent  the  applied  pressures  where  P^,  Pm> 
and  Pp  are  the  external  pressures  applied  to  the  triangular  faces  opposite 
to  nodes  j,  m,  and  p respectively.  The  pressures  are  all  multiplied  by  a 


common  time-dependent  factor,  Ffc.  The  terms  in  the  second  sets  of  paren 
theses  represent  the  effect  of  the  internal  element  stresses.  The  forces  at 


the  other  nodes  are  readily  obtained  by  changing  subscripts.  The  net  forces 

at  node  i (F1,  f'  , F1 ) are  the  sum  of  the  node  i forces  from  each  element 
x y z 

which  includes  that  node. 


6.  EQUATIONS  OF  MOTION 

The  equations  of  motion  can  be  numerically  integrated  by  assuming  a con- 
stant velocity  for  each  time  increment.  The  acceleration  of  node  1 in  the 
x direction  at  time  = t is 

*\-ir  <SI» 

i 

The  new  constant  velocity  for  the  next  time  increment  is 

ut+  = m + u * At  (52) 

i i 1 

where  u|”  is  the  constant  velocity  for  the  previous  time  increment  and  At 
is  the  average  of  the  two  integration  time  increments  about  time  = t. 
Finally,  the  new  displacement  at  time  = t + At  is 

t + At  t • t+  , / c o \ 

u.  = u-  + u.  At  too; 

1 L l 

The  equations  of  motion  for  the  y and  z directions  have  a similar  form. 
The  integration  time  increment  used  for  the  equations  of  motion  is  given 


At  = C, 


/ 2 ^ 2 

(g  + CQ 


where  - CUQ/p,  h is  the  previously  determined  minimum  altitude,  and 
Q 

c is  the  sound  velocity.  The  constant,  C,,  must  be  less  than  unity  to  en- 
s 

sure  that  At  is  always  less  than  the  time  required  to  travel  across  the 
shortest  dimension  of  the  tetrahedron  at  the  sound  velocity  of  the  material. 
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Use  of  a time  increment  significantly  larger  than  that  specified  by  Equation 
(54)  will  lead  to  numerical  instability  (Reference  5).  The  EPIC-3  program 
restricts  the  time  increment  from  increasing  more  than  10  percent  per 
cycle. 

7.  SLIDING  SURFACES 

It  is  sometimes  necessary  to  allow  for  sliding  to  occur  between  two  sur- 
faces. A summary  of  the  important  steps  for  the  sliding-surface  technique 
follows : 

• Identify  a master  sliding  surface  defined  by  an  orderly 
arrangement  of  master  nodes. 

• Identify  slave  nodes  which  may  slide  along  the  master 
surface. 

• For  each  integration  time  increment,  apply  the  equations 
of  motion  to  both  the  master  nodes  and  the  slave  nodes  in 
the  usual  manner. 

• For  each  slave  node,  find  the  triangular  plane  (defined  by 
three  master  nodes)  whose  projection  contains  that  slave 
node. 

• Check  to  determine  if  there  is  interference  between  the 
slave  node  and  the  master  surface  (the  triangular  plane 
defined  by  three  master  nodes). 

• If  there  is  interference,  place  the  slave  node  on  the  master 
surface  in  a direction  normal  to  the  appropriate  master 
triangular  plane. 
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• Calculate  the  momentum  change  (impulse)  to  the  slave  node 
caused  by  placing  it  on  the  master  surface;  then  transfer 
this  lost  momentum  to  the  three  surrounding  master  nodes 
in  a direction  normal  to  the  master  surface. 

• Calculate  the  frictional  impulse  developed  between  the  slave 
node  and  the  three  master  nodes.  The  impulse  is  applied  in 
the  plane  of  the  master  surface  and  in  a direction  opposing  the 
relative  motion. 

The  first  step  in  the  process  is  to  define  the  master  surface.  Jt  is  generally 
desirable  to  define  this  surface  with  an  orderly  arrangement  of  nodes  since 
a systematic  search  is  required  on  this  surface.  An  example  of  this  is 
shown  at  upper  left  in  Figure  3 where  a flat  plate  is  represented  by  a finite- 
element  model  of  nodes  and  elements.  The  z axis  points  vertically  upward 
and  the  x and  y axes  are  in  a horizontal  plane.  Although  the  basic  principles 
involved  in  this  technique  are  independent  of  the  orientation  of  the  axes,  the 
search  technique  is  simplified  if  the  z coordinates  of  the  master  surface  are 
limited  to  a single-valued  function  of  x and  y (i.  e.  , any  vertical  line  par- 
allel to  the  z axis  must  not  pass  through  the  master  surface  at  more  than 
one  point).  This  specialized  case  will  be  described  first,  and  the  specific 
characteristics  of  the  generalized  case  will  be  described  later. 

P’or  the  specialized  case,  the  slave  surface  is  above  the  master  surface 
relative  to  the  z axis.  Due  to  the  symmetry  of  the  plate  about  a plane  con- 
taining the  x and  z axes,  only  one-half  the  plate  is  considered  as  shown.  The 
EPIC-3  code  uses  a master  surface  v.efined  by  M rows  of  nodes,  with  each 
row  containing  N nodes.  Looking  downward  from  the  positive  z axis,  the 
second  row  of  nodes  from  N + 1 to  2N  is  to  the  left  of  the  first  row  of  nodes 
going  from  1 to  N.  Likewise,  the  third  row  is  to  the  left  of  the  second  row, 
etc.  For  the  specific  case  shown  in  Figure  3,  N = 17  and  M = 4,  for  a total 
of  68  master  nodes. 
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Figure  3.  Sliding  Surface  Procedure 
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It  is  also  necessary  to  identify  slave  nodes  for  the  other  sliding  surface. 
Since  these  nodes  do  not  require  an  orderly  arrangement,  it  is  usually  con- 
venient to  designate  the  more  geometrically  complex  surface  as  the  slave 
surface.  It  is  also  desirable  to  restrict  the  slave  node  spacing  from  becom- 
ing significantly  greater  than  the  master  node  spacing,  because  this  could 
introduce  localized  deformations  in  the  master  surface  at  the  slave  node 
locations.  Likewise,  the  mass  of  the  slave  nodes  should  not  be  significantly 
greater  than  the  mass  of  the  master  nodes,  since  this  could  result  in  un- 
realistically high  velocities  of  the  master  nodes  at  impact. 

For  each  integration  loop,  the  equations  of  motion  are  applied  to  both  the 
master  and  slave  nodes  in  the  manner  previously  described.  It  is  then  nec- 
essary to  check  each  slave  node  to  determine  if  it  has  passed  through  the 
master  surface,  thus  causing  interference.  To  begin,  the  extreme  values 
of  the  entire  master  surface  are  determined;  if  the  slave  node  is  not  within 
this  region,  there  can  be  no  interference.  If,  for  instance,  the  z coordinate 
of  the  slave  node  is  greater  than  the  maximum  z coordinate  of  the  master 
surface,  there  can  be  no  interference,  and  the  check  for  that  particular 
slave  node  is  complete. 

If  the  slave  node  is  within  the  confines  of  the  master  surface,  a search 
through  the  master  surface  must  be  initiated.  It  is  for  this  search  that  it 
is  desirable  to  have  an  orderly  arrangement  of  master  nocfes.  The  lower 
left  portion  of  Figure  3 shows  the  projection  of  an  arbitrarily  displaced 
master  surface  on  the  x-y  plane.  The  triangular  grid  is  consistent  with  the 
triangular  faces  of  the  tetrahedron  elements.  Since  z is  a single- valued 
function,  there  is  no  overlapping  of  triangles  on  the  x-y  projection.  The 
objective  of  the  search  is  to  determine  which  triangle  contains  the  x-y  pro- 
jection for  each  slave  node.  The  search  begins  by  considering  the  tri- 
angles between  the  first  two  rows  of  nodes.  With  reference  to  the  lower 
left  portion  of  Figure  3,  the  first  triangle  is  defined  by  nodes  1,  2,  18,  the 
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second  by  nodes  2,  19,  18,  and  the  third  by  nodes  2,  3,  19,  etc.  If  the 
slave  node  projection  is  not  found  by  the  last  triangle  in  the  row  (nodes  17, 
34,  33),  the  pattern  is  repeated  in  the  second  row  of  triangles  beginning 
with  the  triangle  defined  by  nodes  18,  19,  35, 


An  efficient  way  to  check  if  the  slave  node  projection  is  within  a specific 
triangle  is  to  determine  the  distance  from  the  slave  node  to  each  of  the 
lines  defining  the  three  sides  of  the  triangle.  The  distance  to  line  i-j  in 
the  x-y  plane  is  given  by 

£lx's  +32y/s  +/33 


xy 


(55) 


where  /3  = y.  - y (3  = x.  - x.  and  = x y . - x y.,  The  coordinates  of 

1 1 J Z J 1 ° 1 J J 1 

master  node  i are  x^,  y^  and  z^,  and  the  coordinates  of  the  slave  node  are 

x/g»  y,g>  and  z>s-  ^ master  node  j is  counterclockwise  from  node  i,  when 

looking  downward  from  the  positive  z axis,  and  6 is  positive  for  all  three 

xy 

lines  of  the  triangle,  then  the  slave  node  is  contained  within  the  triangle. 
This  is  also  illustrated  in  Figure  3 where  the  slave  node  is  included  in  the 
triangle  defined  by  nodes  31,  32  and  48. 


After  the  appropriate  triangle  is  identified,  it  is  necessary  to  determine  if 
there  is  interference.  The  equation  of  a plane  through  nodes  i,  j and  k 
(taken  counterclockwise,  when  looking  downward  from  the  positive  z axis) 
is 


Ax  + By  + Cz  + D = 0 


(56) 


where 
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(57a) 


(57  b) 


(57c) 


(57d) 


In  the  upper  right  portion  of  Figure  3,  a vector  is  shown  which  passes 
through  the  initial  slave  node  position  and  is  normal  to  the  triangular  plane. 
The  three  direction  cosines  of  this  vector  are  readily  obtained  by  dividing 
A,  B and  C by  \Ja^  + B^  + C^.  Subsequent  formulations  will  consider  the 
constants  (A,  B,  C,  D)  to  be  redefined  such  that  A,  B,  C represent  the 
direction  cosines.  The  normal  distance  between  the  slave  node  and  the 
master  plane  is 


S = ' (Axs  + By's  + Cz's  + D)  (58) 

where  x'  y ' and  z ' are  the  coordinates  of  the  slave  node.  If  6^T  is 

S S S IN 

positive,  there  is  interference  and  the  slave  node  is  placed  onto  the  master 
surface  in  a direction  normal  to  the  master  surface.  The  resulting  final 
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coordinates  of  the  slave  node  are  x = x'  + Ax,  y = y'  + Av  and  z - z' 

s s ’ J s J s 7 s s 

+ Az,  where  Ax,  Ay  and  Az  are  obtained  by  multiplying  5^.  by  the  appropriate 
direction  cosines  of  the  normal  vector. 

After  geometric  compatibility  has  been  achieved  by  placing  the  slave  node 
on  the  master  surface,  it  is  necessary  to  adjust  the  equations  of  motion  to 
account  for  this  change.  Since  the  geometric  changes  are  made  only  in  the 
direction  normal  to  the  master  plane,  the  affected  equations  of  motion  are 
altered  only  in  this  normal  direction.  It  should  be  noted  that  frictional 
effects  alter  the  equations  of  motion  in  the  plane  of  the  master  surface,  but 
these  will  be  considered  later. 

In  the  process  of  moving  the  slave  node  to  the  master  surface,  an  effective 
velocity  change  is  imposed  on  the  slave  node.  This  loss  of  velocity  in  the 
normal  direction  is  AV^  = 6^/  At,  where  At  is  the  integration  time  incre- 
ment for  the  previous  cycle.  This  velocity  change  also  gives  a momentum 
change  along  the  normal  direction  equal  to  M A\\,,  where  M is  the  mass 
of  the  slave  node.  This  momentum  is  then  transferred  to  the  three  surround- 
ing master  nodes.  If  R.,  R.  and  R,  represent  the  fractions  of  the  momentum 

1 3 K 

to  be  transferred  to  the  three  nodes,  a summation  of  the  translational  and 
rotational  momenta  of  the  three  nodes  gives  the  following  equations: 

Rl  + R + Rk  = 1 (59) 

Ri  <xi  ' V + Rj  (xj  ' xs>  + Rk  (xk  “ xs)  = 0 (60) 

Ri  - ys)  + Rj  (yj  - ys)  + Rk  <yk  - = 0 {61) 

These  three  equations  and  three  unknowns  define  the  appropriate  momen- 
tum distribution  to  the  three  nodes.  This  momentum  is  applied  to  the 
master  nodes  in  the  form  of  instantaneous  velocity  changes,  parallel  to  the 
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normal  vector  but  in  the  opposite  direction.  For  instance,  the  normal 
direction  velocity  change  to  node  i is 

AVi  = RiMsAVN/Mi  (62) 

where  NL  is  the  mass  of  master  node  i and  Mg  is  the  mass  of  the  slave 
node.  As  was  done  for  the  displacements,  the  specific  velocity  changes 
in  the  x,  y and  z directions  are  obtained  by  multiplying  AV-  by  the  appro- 
priate direction  cosines. 


It  is  sometimes  desirable  to  include  frictional  effects  in  the  analysis.  This 
effect  is  constrained  to  the  plane  of  the  master  surface  and  the  direction 
opposes  the  relative  motion.  The  net  magnitude  of  the  frictional  velocity 
change  to  the  slave  node  is  proportional  to  the  normal  velocity  change  to 
the  slave  node,  or 


iVP  ’ VVN 


where  f is  the  coefficient  of  friction. 


(63) 


To  apply  this  velocity  change  in  the  direction  of  the  relative  velocities  it 
is  necessary  to  determine  the  velocity  of  the  master  surface  at  the  slave 
node  location.  The  procedure  is  similar  to  that  used  in  Equations  (1) 
through  (4)  which  determined  the  velocity  distribution  within  a tetrahedron. 
For  this  case,  however,  the  velocity  is  determined  within  a two-dimensional 
triangle  since  it  can  be  defined  as  a function  of  the  x and  y coordinates. 

The  resulting  master  surface  velocity  in  the  x direction,  at  the  slave  node 
position,  is 


Vi  + Vj 


+ qkuk 


(64) 
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The  other  velocity  components  (v  and  wm)  have  a similar  form  to  Equation 


(64),  and  the  geometry  constants  (q-,  q.,  q.  ),  have  the  form 

1 3 K 

qi = 2 x“  cx.yk  - Vj + (yj  - yk>  xs  + (xk  - V ys]  (65) 


xy 


where  A r is  the  projected  cross-sectional  area  of  the  triangular  master 

xy 


surface  on  the  x-y  plane,  and  the  other  terms  represent  the  coordinates  of 
the  master  and  slave  nodes. 


Now  that  the  velocities  of  both  the  master  surface  and  slave  node  are  defined, 
it  is  necessary  to  determine  the  components  of  these  velocities  which  are  in 
the  plane  of  the  master  surface.  This  is  accomplished  by  subtracting  the 
velocity  components  normal  to  the  master  plane,  from  the  total  velocity 
components.  These  velocity  components  for  the  slave  node  are  expressed 
as 


c • 

II 

c • 

-AVI? 

(66) 

s s 

s 

, p . 

V = V 

N 

- B VG 

(67) 

s s 

s 

, P 

w„  = w„ 
s s 

(68) 

where  the  slave  node  velocity  normal  to  the  master  surface  is 


,N 


Va  = Au  + Bv  + Cw 
s s s s 


(69) 


. p • p . p 

The  in-plane  velocities  of  the  master  surface  (u  , vm,  wm)  are  obtained 


in  a similar  manner. 


J 
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The  relative  velocity  components  (u^,  vj^,  are  then  obtained  in 

the  form  of 


• P _ * P • P 

rel  Um  Us  (70) 


I'  inally,  the  friction  induced  velocity  changes  (Au^*, 
form 


Av 


Aw  ) have  the 


rel 


AV, 


f(urel)2  +(;rel)2  + (wreir 


(71) 


The  momentum  change  of  the  master  nodes,  due  to  the  frictional  force, 
gives  m-plane  velocity  changes  to  master  node  i (Au^  Av^  Aw^ ) of  the  form 

Pi  \/M.\ 

(72) 


This  completes  the  description  of  the  sliding  surface  procedure  for  the 
specialized  search  technique  shown  in  Figure  3.  The  necessity  for  a gen- 
eralized search  technique  is  shown  in  the  lower  right  portion  of  the  same 
figure  where  the  master  surface  is  not  a single-valued  function.  A slave 
node  is  shown  slightly  below  the  master  triangle  containing  nodes  7 and  8, 
and  it  should  properly  be  placed  on  this  triangle.  If  the  specialized  search 
technique  would  be  used  for  this  condition,  however,  the  slave  node  would 
be  placed  on  the  triangle  containing  nodes  2 and  3 since  it  is  also  contained 
in  the  x-y  projection  of  this  triangle. 


To  alleviate  this  problem,  the  generalized  search  technique  selects  the 
master  surface  triangle  which  is  closest  to  the  slave  node.  This  selected 
:r  .angle  must  also  be  within  a zone  of  acceptability  which  is  dependent  on 


28 


the  dimensions  of  the  selected  triangular  surface.  If  it  is  within  this  zone 
a check  for  interference  is  made  with  Equation  (58).  Placement  of  the  slave 
node  onto  the  master  surface  and  the  adjustments  of  the  nodal  velocities 
are  done  in  a similar  manner  to  that  which  was  done  for  the  specialized 
formulation.  In  some  cases,  however,  the  master  triangle  must  be  defined 
in  terms  of  the  x-z  or  y-z  coordinate  system,  depending  on  the  orientation 
of  the  triangle. 

The  specialized  search  routine  is  computationally  faster  than  the  general- 
ized search  routine  since  it  generally  does  not  search  through  every  tri- 
angle on  the  master  surface;  if  the  slave  node  is  beyond  the  extremes  of  the 

master  surface,  the  routine  does  not  even  begin  the  search,  and  if  it  does 
search,  it  stops  when  the  proper  triangle  is  located.  For  either  routine  the 
translational  momenta  are  absolutely  conserved  in  all  three  directions. 

Small  errors  are  introduced  into  the  center  of  gravity  positions  when  the 
slave  node  is  moved  to  the  master  surface,  since  there  is  no  corresponding 
movement  of  the  master  nodes,  which  receive  only  instantaneous  velocity 
changes.  This  center  of  gravity  error  also  introduces  small  errors  into 
the  rotational  momenta. 
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SECTION  III 

COMPUTER  PROGRAM  DESCRIPTION 


1.  BACKGROUND 

This  section  describes  the  important  characteristics  of  the  computer  pro- 
gram. The  EPIC-3  software  was  developed  in  FORTRAN  on  a Honeywell 
series  6000  large-scale  digital  computer.  The  program  is  quite  portable; 
only  minor  changes  were  required  to  make  it  operational  on  the  Eglin 
Air  Force  Base  CDC  6600  computer.  Two  major  design  goals  were  to  de- 
velop the  program  for  efficient  execution  on  fourth-generation  vector  com- 
puters and  to  have  the  capability  to  provide  solutions  for  problems  of  un- 
limited size. 

Even  though  the  advanced  EPIC -3  code  is  a descendant  of  an  earlier  version, 
it  was  totally  redesigned  and  rewritten.  The  top-down,  hierarchical  design 
approach  resulted  in  a structured  collection  of  functional  modules.  The  soft- 
ware was  implemented  with  an  emphasis  on  clarity  and  readability.  With 
few  exceptions,  strict  coding  conventions  and  a limited  set  of  coding  con- 
structs were  employed.  The  design  documentation,  written  as  comment  lines, 
has  been  embedded  into  the  source  code  for  each  module.  Vectorization 
guidelines  were  applied  to  major  computational  loops,  whenever  possible,  to 
enhance  execution  speed  on  computers  that  have  vector-processing  capability. 
Extensive  use  of  files  for  intermediate  storage  was  required  to  provide  solu- 
tions for  large  problems  which  cannot  be  core  contained. 

The  EPIC-3  Postprocessor  provides  the  user  with  a flexible  plot  capability 
to  display  two-  or  three-dimensional  views  of  the  initial  and  deformed  geo- 
metry as  well  as  plots  of  velocity  vectors.  The  three-dimensional  capa- 
bility is  significant  because  any  perspective  is  available,  and  hidden  lines 
can  be  optionally  removed  in  the  geometry  plots. 
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2.  PROGRAM  STRUCTURE 


The  EPIC -3  software  is  currently  implemented  on  an  H6080  computer.  It 
utilizes  the  H6080  system  sort  package  and  standard  CALCOMP  Plot  calls 
in  addition  to  intrinsic  functions  supplied  with  its  FORTRAN  complier.  The 
software  is  also  operational  on  the  Eglin  AFB  CDC  6600  computer,  which 
uses  a Stromberg-Carlson  Plotter. 

As  shown  in  Figure  4,  EPIC-3  has  three  separate  programs  which  communi- 
cate by  means  of  one  restart  file  format.  The  Preprocessor  (PREP)  sets  up 
the  initial  problem  data,  writing  it  on  the  restart  file.  The  Main  Routine 
(MAIN)  reads  this  restart  file  or  one  created  by  a previous  Main  Routine  run, 
and  then  it  processes  additional  time  cycles,  outputting  to  a new  restart  file 
at  user-specified  intervals.  The  Postprocessor  (POST)  obtains  its  data  from 
restart  files  created  by  either  the  Main  Routine  or  the  Preprocessor. 

The  three  programs  were  functionally  decomposed  with  a top-down  approach 
to  design.  The  resultant  program  structures  and  module  functions  are  pre- 
sented as  Figures  5 through  7.  The  data  handling,  which  is  common  to  all 
EPIC -3  programs,  is  facilitated  by  the  multi-use  service  routines  listed  in 
Figure  8.  A summary  of  all  EPIC-3  subroutines,  with  their  associated  com- 
mon blocks  and  service  routines,  is  provided  in  Figure  9. 


3.  VECTORIZATION  GUIDELINES 

A requirement  for  EPIC -3  is  that  it  will  eventually  be  used  on  vector  com- 
puters such  as  the  CRAY,  and  it  must  therefore  be  written  to  be  compatible 
with  these  machines.  Compilers  on  vector  computers  such  as  the  CRAY 
analyze  the  innermost  loops  of  a program  and  attempt  to  use  vector  instruc- 
tions which  exploit  the  advanced  hardware  and  software  features  of  the  com- 
puter. Since  no  special  forms  are  required  for  vector  operations,  the 
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Figure  5.  Preprocessor  Hierarchy  Chart 
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Figure  6.  Main  Routine  Hierarchy  Chart 
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Figure  7. 


Postprocessor  Hierarchy  Chart 
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SERVICE 
ROUTI  NE 


DESCRIPTION 


EREAD  READS  ELEMENT  DATA  FROM  FILE  TO  CORE 

EREADI  READS  ELEMENT  DATA  FROM  FILE  TO  BUFFER 

EREAD2  READS  ELEMENT  DATA  FROM  BUFFER  TO  CORE 

ERTAPE  READS  ELEMENT  DATA  FROM  RESTART  TAPE  TO  CORE 

EWRITE  WRITES  ELEMENT  DATA  ONTO  FILE  FROM  CORE 

EWTAPE  WRITES  ELEMENT  DATA  ONTO  RESTART  TAPE  FROM  CORE 

HDATA  PRINTS  INPUT  DATA 

LBPOS  DETERMINES  CORE  POSITIONS  OF  ELEMENT  BLOCK 
LPOS  DETERMINES  CORE  POSITION  OF  INDIVIDUAL  ELEMENT 

NBPOS  DETERMINES  CORE  POSITIONS  OF  NODE  BLOCK 
NFIX  DECODES  IXYZ  ARRAY  FOR  NODE  INDICATORS 

NPOS  DETERMINES  CORE  POSITION  OF  INDIVIDUAL  NODE 

NREAD  READS  NODE  DATA  FROM  FILE  TO  CORE 

NREADl  READS  NODE  DATA  FROM  FILE  TO  BUFFER 

NREAD2  READS  NODE  DATA  FROM  BUFFER  TO  CORE 

NRTAPE  READS  NODE  DATA  FROM  RESTART  TAPE  TO  CORE 

NWRITE  WRITES  NODE  DATA  ONTO  FILE  FROM  CORE 

NWTAPE  WRITES  NODE  DATA  ONTO  RESTART  TAPE  FROM  CORE 

PERSPC  TRANSFORMS  3-0  COORDINATES  TO  2-0  PERSPECTIVE 
SWITCH  REWINDS  AND  INTERCHANGES  TWO  FILES 

Figure  8.  Service  Routines 
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Figure  9.  EPIC-3  Subroutines 
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NRTAPE 


NWTAPE 


transportability  of  a vector -oriented  FORTRAN  program  to  other  compilers 
is  not  endangered.  During  execution  on  a vector  computer,  a vector  loop 
can  run  many  times  as  fast  as  its  equivalent  scalar  loop.  However,  there 
are  certain  problems  that  prevent  a loop  from  vectorizing. 

To  be  vectorizable , a computational  loop  should  manipulate  arrays  and  store 
the  results  of  computations  in  arrays.  Loops  that  contain  branches  such  as 
GO  TO's  and  IF's,  cannot  currently  be  vectorized.  The  use  of  temporary- 
variables,  CALL  statements  or  pointers  (an  array  position  identifier  ob- 
tained from  another  array)  also  prohibits  vectorization. 

The  following  general  rules  apply  for  vector  loops:  , 

• Keep  subscripts  simple  and  explicit. 

• Do  not  use  temporary  variables,  pointers,  IF's,  GO  TO's, 
or  CALL'S. 

• Use  mathematical  functions  for  which  there  are  vectorized 
versions  in  the  target  vector  computer. 

• If  a large  loop  contains  only  a few  unvectorizable  statements, 
try  to  rewrite  it  as  two  or  more  loops,  with  at  least  one  being 
vectorizable . 

Much  of  the  EPIC-3  code  is  currently  vectorized  based  on  the  application  of 
the  aforementioned  rules.  Since  most  of  the  element  computations  depend 
on  nodal  data  such  as  positions  and  velocities,  the  nodal  data  are  redefined 
in  terms  of  element  variables  such  that  pointers  are  not  required  for  the 
primary-element  computational  loops.  The  gathering  of  nodal  data  into  ele- 
ment arrays  is  performed  in  subroutine  EGET,  and  the  subsequent  distribu- 
tion of  element  data  back  to  the  nodal  arrays  (concentrated  forces)  is  per- 
formed in  subroutine  EPUT. 
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4.  CODING  GUIDELINES 

Most  of  the  EPIC -3  program  (with  the  exception  of  the  geometry  generators) 
is  written  using  structured  programming  guidelines.  Other  exceptions  occur 
in  some  of  the  primary  computational  loops  where  computational  efficiency 
or  vectorization  would  be  impaired  if  the  structured  programming  guidelines 
were  incorporated.  Every  attempt  has  been  made  to  provide  consistency  and 
clarity  in  all  the  subroutines.  An  example  of  the  coding  is  provided  with  sub- 
routine MOTION,  shown  in  Figure  10. 

5.  DATA  STRUCTURE 

The  significant  data  within  the  EPIC -3  software  are  organized  in  labeled  com- 
mon blocks.  These  blocks  provide  the  basic  means  for  communication  be- 
tween subroutines.  The  block  labels,  in  addition  to  the  variable  names 
themselves,  are  intended  to  be  as  self-explanatory  as  possible.  Figure  9 
lists  the  common  block  usage  by  subroutine.  The  common  blocks  and  their 
contents  are: 

/BOOKKP/  - All  bookkeeping  information  such  as  block  sizes,  num- 
ber of  blocks,  total  number  of  elements,  and  total  num- 
ber of  nodes . 

/EBUFR/  - Buffer  for  element  data  that  are  being  read. 

/EBUFW/  - Buffer  for  element  data  that  will  subsequently  be  written 

onto  an  element  file. 

/ELEMNT/  - Data  arrays  for  one  block  of  elements. 

/FILES/  - Variables  used  to  store  the  logical  unit  designations  for 
all  EPIC -3  files. 
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Figure  10.  Example  of  Structured  Programming  Coding  Format 


/MATERL/ 
/MISC/ 
/NBUFR / 
/NBUFW  / 

/NODE/ 

/PLOTD/ 

/SLIDED/ 

/VECTOR/ 


Material  properties. 

Miscellaneous  data  needed  by  multiple  routines. 

Buffer  for  node  block  data  that  are  being  read. 

Buffer  for  node  block  data  that  will  subsequently  be  writ- 
ten onto  a node  file. 

Data  arrays  for  the  entire  bandwidth  of  node  blocks. 
Assorted  plot  data. 

Sliding  surface  data. 

Arrays  needed  for  various  element  vector  loops. 


Files  are  used  extensively  by  the  EPIC-3  software.  The  node,  element,  and 
sliding  surface  files  are  designed  to  minimize  the  program's  core  storage  re- 
quirements by  providing  for  the  storage  and  subsequent  retrieval  of  data  that 
are  not  currently  needed  for  calculations.  Typically,  data  will  be  read  from 
files  into  common  block  locations,  it  will  be  processed,  and  then  the  updated 
data  will  be  written  to  other  files.  A summary  of  the  various  files  is  given 
in  Section  IV. 


6.  BANDWIDTH  DETERMINATION 

Since  three-dimensional  problems  often  require  many  nodes  and  elements, 
it  is  not  generally  possible  to  contain  all  data  within  the  central  memory  core 
of  existing  computers.  As  a result,  it  is  necessary  to  use  disk  files  for  in- 
termediate storage.  Since  there  is  interaction  between  various  nodes  for 
sliding  surface  computations,  and  since  element  computations  require  various 
node  data,  it  is  necessary  to  provide  a logical,  orderly  technique  whereby  the 
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appropriate  element  and  node  data  are  in  core  at  the  proper  time  during  the 
computations . 

For  each  cycle  of  numerical  integration,  there  are  three  computational  loops. 
A brief  discussion  of  each  follows: 

• The  first  loop  computes  the  equations  of  motion  of  the  nodes 
(Subroutine  MOTION).  It  does  not  include  the  sliding  sur- 
face computations.  The  computations  for  each  node  are  self 
contained;  they  do  not  depend  on  data  from  other  nodes  or 
elements.  Therefore,  this  computation  has  no  effect  on  the 
data  handling  or  central  memory  requirements. 

• The  second  loop  adjusts  the  equations  of  motion  for  the  slid- 
ing surface  nodes  (Subroutine  SLIDE).  The  computations 
for  each  slave  node  do  not  depend  on  any  element  data;  they 
do,  however,  interact  with  data  from  other  nodes  (the  mas- 
ter nodes).  Therefore,  for  a specific  sliding  surface,  it  is 
necessary  that  all  master  nodes  be  in  core  when  each  slave 
node  is  processed. 

• The  third  loop  computes  various  data  for  the  elements  (Sub- 
routine ELOOP).  The  computations  for  each  element  do  not 
depend  on  data  from  any  other  element.  They  do,  however, 
depend  on  node  data.  The  computations  for  a specific  ele- 
ment require  the  position  and  velocity  data  for  the  four  nodes 
which  define  that  element.  The  element  then  uses  this  data 
(together  with  material  property  data)  to  generate  concen- 
trated forces  which  it  distributes  to  the  four  nodes.  There- 
fore, for  a specific  element  being  processed,  it  is  neces- 
sary for  the  four  nodes  of  an  element  to  be  in  core  as  that 
element  is  being  processed. 
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The  second  and  third  loops  both  place  requirements  on  the  data  handling  and 
central  memory  requirements.  It  is  from  these  requirements  that  the  con- 
cept of  a bandwidth  is  introduced.  The  third  loop  is  the  most  complicated 
since  both  node  and  element  data  must  be  in  central  memory  at  the  proper- 
time.  The  discussion  which  follows  is  for  this  third  (element)  loop. 

An  unordered  arrangement  of  nodes  and  elements  is  shown  at  upper  left  in 
Figure  11.  Each  element  along  the  vertical  axis  is  affected  by  four  nodes 
along  the  horizontal  axis.  If  the  primary  computational  activity  is  based  on 
processing  element  data  in  a sequential  manner,  the  unordered  arrangement 
does  not  restrict  the  location  of  the  appropriate  four  nodes  for  that  element. 
As  a result,  if  the  required  node  data  are  not  in  core  but  must  be  brought  in 
from  disk  storage  in  a random  manner,  the  computations  are  inefficiently 
performed. 

In  contrast,  an  ordered  arrangement  of  nodes  and  elements  is  shown  at  upper- 
right  in  Figure  11.  Here  it  can  be  seen  that  the  four  nodes  for  each  element 
are  restricted  to  fall  within  a range  defined  as  the  bandwidth.  This  means 
that,  as  each  element  is  processed,  it  is  necessary  to  provide  access  to  only 
a limited  number  of  nodes  contained  within  the  bandwidth.  This  ordered  ar- 
rangement is  obtained  by  numbering  the  nodes  and  elements  in  a logical  and 
consistent  manner.  Generally,  the  bandwidth  for  the  element  loop  is  mini- 
mized if  the  nodes  are  increasingly  numbered  in  a direction  normal  to  the 
cross  section  with  the  smallest  number  of  nodes. 

The  lower  portion  of  Figure  11  shows  an  ordered  arrangement  of  nodes  and 
elements  where  the  data  are  represented  by  N blocks  of  node  data  and  M 
blocks  of  element  data.  These  blocks  contain  data  for  groups  of  nodes  and 
elements.  This  is  done  to  simplify  the  transfer  of  data  between  disk  files 
and  central  memory  and  also  to  allow  the  vectorized  computational  loops  to 
act  on  an  entire  block  of  data.  The  bandwidth  for  a particular  element  block 
is  the  number  of  consecutive  node  blocks  needed  to  contain  both  the  lowest 
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Figure  11.  Bandwidth  Representation  for  Node- 
Element  Connectivity 


nude  number  and  the  highest  node  number  used  by  any  element  in  that  element 
block  or  any  preceding  block.  A bandwidth  of  five  node  blocks  is  shown  to 
■ ontain  the  require  I node  data  for  each  element  block  in  Figure  11.  The  node 
i.jck>  represented  with  solid  lines  are  those  blocks  which  are  needed  by  the 
>r responding  element  block,  and  the  blocks  represented  with  dashed  lines 
u . nut  nee  le  I by  the  element  blocks.  It  can  be  seen  that  this  bandwidth  is 
letermine  j b\  element  block  J,  since  it  requires  access  to  the  five  node 
tdocks  iron  block  1-2  to  block  l*-2. 

\nuther  requirement  for  the  bandwidth  comes  from  the  sliding  surface  compu- 
tation.-. Since  this  involves  only  the  nodes,  it  is  not  as  complicated  as  the 
previously  defined  bandwidth  obtained  from  the  elements.  For  a specified 
sliding  surtnce,  it  the  slave  nodes  all  have  lower  node  numbers  than  the  as- 
sociated  nodes  of  the  nia-ter  surface,  then  the  bandwidth  is  the  number  of 
consecutive  node  blocks  needed  to  contain  ail  slave  and  master  nodes  for  that 
particular  sliding  surface.  The  final  required  bandwidth  is  the  maximum  de- 
termined from  the  elements  and  the  sliding  surfaces.  Since  it  is  desirable 
to  minimize  this  bandwidth,  the  geometry  generators  in  the  Preprocessor 
are  developed  aci  or  iingly. 

1 he  current  vti  sion  >1  t'.i  code  contain.-'  64  nodes  per  node  block  and  G4  ele- 
ments per  element  block.  (The  node  and  element  block  sizes  do  not  have  to 
be  identical,  however.  ) 1'his  size  is  well  suited  for  vector  processing  on  the 
CliA'i  computer,  li  also  provides  fot  the  transfer  of  a significant  amount  of 
data  such  that  the  tune  required  to  initiate  the  Input/Output  activity  is  rela- 
tively small.  It  should  be  noted  that  some  problems  will  have  few  enough 
nodes  so  all  node  data  can  be  core-contained.  If  adequate  core  is  provided, 
the  transfer  of  nodal  data  between  files  and  core  will  be  eliminated.  The 
savings  in  the  associated  Input/Output  and  Central  Processor  time  for  the 
core-contained  option  must  be  compared  to  the  advantages  associated  with 
using  a smaller  core  size  which  contains  only  a bandwidth  of  nodal  data. 
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7.  COMPUTATIONAL  SEQUENCE 

The  primary  computational  sequences  for  EPIC- 3 have  been  designed  to 
minimize  problem  size  requirements  and  to  enhance  vector  processing.  A 
schematic  representation  of  the  data  manipulation  is  shown  in  Figure  12. 

The  central  memory  core  contains  one  block  of  element  data  and  a bandwidth 
of  node  blocks.  Node  and  element  buffers  for  read  and  write  activity  are  re- 
quired, and  storage  space  is  provided  for  the  material  and  sliding  surface 
data. 
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Figure  12.  Program  Logic  and  Memory  Requirements 


The  following  process  describes  the  important  steps  which  occur  during  an 
element  loop  (Subroutine  ELOOP).  These  steps  are  illustrated  in  Figures 
11  and  12  where  the  bandwidth  is  five  node  blocks. 


• Initial  Step  to  Set  Up  Processing  - A complete  bandwidth  of 
node  blocks  is  read  into  core  from  the  current  node  read 
file.  An  additional  node  block,  if  it  exists,  will  be  read  in- 
to the  node  read  buffer.  The  first  element  block  must  also 
be  read  into  core  from  the  element  read  file  and  the  buffered 
read  of  the  next  one  is  initiated.  In  Figures  11  and  12,  this 
would  mean  node  blocks  1 to  5 are  read  into  core,  and  node 
block  6 is  read  into  the  node  read  buffer.  Likewise,  ele- 
ment block  1 is  read  into  core,  and  element  block  2 is  read 
into  the  element  read  buffer. 

• Basic  Processing  Sequence  (Done  for  Each  Element  Block)  - 
The  following  must  be  repeatedly  done  until  the  maximum 
node  block  needed  for  this  element  block  is  in  core: 

The  lowest  numbered  node  block  in  core  is  trans- 
ferred to  the  node  write  buffer  and  written  on  the 
node  write  file.  Since  it  is  no  longer  within  the 
bandwidth,  it  will  not  be  affected  by  subsequent  ele- 
ments. Then  the  next  node  block  is  transferred  from 
the  node  read  buffer  to  the  vacated  core  space  and 
another  buffered  node  read  is  begun.  In  Figures  11 
and  12,  before  element  block  J can  be  processed, 
node  block  1-3  must  be  replaced  with  node  block  1+2 

The  element  computations  for  this  block  can  now  be  per- 
formed since  the  required  node  data  are  in  core.  These  com- 
putations are  summarized  under  Subroutine  ELOOP  in  Figure 
6. 
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The  updated  element  block  is  written  on  the  element  write  file 
and  the  next  element  block  is  transferred  from  the  element 
read  buffer  to  core.  Then  the  next  buffered  element  read  is 
begun.  In  Figures  11  and  12,  element  block  J is  replaced  by 
element  b lock  J+l . 

• Final  Step  to  Complete  the  Processing  - Any  remaining  node 
blocks  in  core  must  be  written  on  the  node  write  file.  This 
would  be  node  blocks  N-4  to  N as  shown  in  Figures  11  and  12. 

The  logical  unit  designators  for  the  node  read  and  write  files 
are  interchanged  so  that  subsequent  processing  will  read  from 
the  updated  file.  Likewise,  the  element  files  are  interchanged. 

The  two  node  loops  (Subroutines  MOTION  and  SLIDE)  are  less  complex  since 
there  is  no  interaction  with  element  data.  For  the  MOTION  loop  the  node 
blocks  are  read  from  one  node  file,  the  computations  are  performed,  and  the 
updated  results  are  written  onto  the  other  node  file. 

For  the  sliding  surface  computations,  a bandwidth  of  node  blocks  must  always 
be  in  core.  If  the  computations  are  initiated  with  the  slave  nodes  in  the  lowest 
node  block,  then  the  affected  higher-numbered  master  nodes  must  always  be 
in  core.  When  the  last  slave  node  in  a specific  sliding  surface  has  been  pro- 
cessed, the  data  (master  node  numbers,  etc.  ) for  the  next  sliding  surface 
are  read  into  core  from  the  sliding  file.  This  requires  the  slave  nodes  of  a 
subsequent  sliding  surface  to  be  numbered  higher  than  those  of  a preceding 
sliding  surface. 


8,  THREE-DIMENSIONAL  PLOTTING  ROUTINES 

This  subsection  summarizes  the  techniques  used  to  obtain  the  three-dimen- 
sional plots.  For  each  three-dimensional  plot  request,  one  call  is  made  to 
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the  SETORG  subroutine  which  develops  the  transformation  from  3-space  to  a 
2-space  plotting  plane.  Since  the  user  specifies  an  eye  point  and  an  origin 
(the  point  looked  at),  SETORG  can  compute  the  equation  for  the  plotting  or 
viewing  plane  which  is  perpendicular  to  the  line  of  sight  and  contains  the  ori- 
gin. Subsequent  calls  to  the  PERSPC  subroutine  provide  2-space  coordinates 
for  any  point  as  it  is  projected  onto  the  viewing  plane. 

The  x,  y,  and  z maximum  and  minimum  values  are  used  to  scale  the  3-D 
plot.  The  eight  combinations  of  coordinates  define  the  corner  nodes  of  a box 
in  space.  PERSPC  is  called  to  project  this  box  onto  the  viewing  plane  using 
the  transformation  defined  for  this  particular  3-D  plot  request.  The  vertical 
dimension  of  the  projected  figure  will  be  scaled  to  fit  within  10  inches  The 
horizontal  dimension  will  use  the  same  scale  to  avoid  distortion;  the  length 
of  the  plot  is  unlimited. 

Plotting  requires  one  complete  pass  through  the  element  file.  All  four  tri- 
angular faces  of  each  element  are  examined.  If  a triangle  fits  within  the 
maximum  and  minimum  extremes,  a corresponding  data  record  is  written 
onto  a new  triangle  file.  The  three  nodes  are  listed  in  ascending  order,  fol- 
lowed by  their  x,  y,  z coordinates.  The  triangle  file  is  sorted  on  ascending 
node  numbers  and  then  a new  file  containing  only  exterior  triangles  is  written. 
The  interior  triangles  are  easily  identified  and  eliminated  by  ttie  presence  of 
two  or  more  consecutive,  identical  records  on  the  sorted  triangle  file. 

The  exterior  triangle  file  is  further  processed  if  the  hidden  line  option  is  re- 
quested. One  complete  node  loop  is  made  to  flag  all  exterior  nodes  (those 
contained  in  the  external  triangles).  To  ascertain  an  exterior  node's  visibility, 
each  exterior  triangle  must  be  examined  to  see  whether  its  2-space  projec- 
tion contains  the  node  and  if  so,  whether  the  node  is  behind  the  triangle. 

Once  all  visible  nodes  are  determined,  one  complete  triangle  pass  is  made 
to  write  those  triangles  having  three  visible  nodes  and  to  create  degenerate 
(two-node)  triangles  when  two  of  the  three  nodes  are  visible.  With  this 
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method,  any  line  of  a triangle  is  eligible  for  eventual  plotting  if  its  endpoints 
are  visible.  Unfortunately,  some  erroneous  lines  can  occur  at  the  corners 
of  a plot  if  the  endpoints  of  a line  are  visible,  but  the  line  itself  is  hidden  by 
another  external  triangle. 

The  final  triangle  file,  with  or  without  hidden  lines,  is  now  processed  for  plot- 
ting. For  each  triangle,  the  PERSPC  subroutine  is  called  three  times  to 
transform  the  three  nodes  to  their  2-space  representation  on  the  plotting 
plane.  Then  standard  plot  subroutines  are  called  to  draw  lines  between  the 
three  nodes . 


SECTION  IV 

PROGRAM  USER  INSTRUCTIONS 


This  section  provides  instructions  for  the  user.  Input  data  for  the  Prepro- 
cessor, the  Main  Routine,  and  the  Postprocessor  are  given  in  the  first  three 
subsections.  Then  the  output  data  are  described  and  instructions  for  chang- 
ing the  dimensions  of  the  program,  the  file  designations,  and  central  proces- 
sor time  estimates  are  given. 

I.  INPUT  DATA  FOR  THE  PREPROCESSOR 

The  functions  of  the  Preprocessor  are  to  define  the  initial  geometry  and  velo- 
city conditions  and  to  determine  the  memory  size  requirements  for  the  Main 
Routine.  A summary  of  the  input  data  for  the  Preprocessor  is  given  in  Fig- 
ure 13,  and  details  of  the  node  and  element  input  data  are  given  in  Figures 
14  and  15.  These  Figures  correspond  to  paragraphs  a,  b,  and  c below.  The 
descriptions  which  follow  are  for  the  data  in  Figure  13.  Consistent  units 
must  be  used. 

a.  Preprocessor  Input  Data  Summary 

• Description  Card  (12A6)  - A description  of  the  problem  pro- 
vided by  the  user. 

• Identification  Card  (315)  - 

CASE  = Case  number  for  run  identification 

0 will  not  print  individual  data  for  each  node  and 
element. 

1 will  print  individual  data 

0 will  not  write  results  on  restart  tape. 

1 will  write  on  restart  tape. 
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Figure  13.  Preprocessor  Input  Data 
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Figure  15.  Element  Input  Data 
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• 22  Material  Cards  (5E15.  8)  - Each  of  the  22  material  cards 

provides  a specific  material  characteristic  for  five  materials. 

For  instance,  the  first  card  describes  the  density  of  the  ma- 
terials, The  density  of  material  1 is  entered  in  columns  1 to  1 5, 
the  density  of  material  2 is  entered  in  columns  16  to  30,  etc.  The 
second  card  contains  the  specific  heats  for  the  five  materials,  etc. 
It  is  only  necessary  to  describe  the  materials  used  for  the 
analysis.  If  materials  1 and  3 are  used  for  a specific  analy- 
sis, all  the  material  1 data  must  be  specified  in  columns  1 to  1 5 
and  all  the  material  3 data  must  be  specified  in  columns  31  to  4 5. 
Columns  16  to  30,  46  to  60,  and  61  to  75,  representing  materials 
2,  4 and  5,  can  be  left  blank. 


Card  1 Density  (mass/ volume) 

Card  2 Specific  heat  (work/mass/degree) 

Card  3 Shear  modulus  of  elasticity  (force/area) 

Card  4 Yield  stress  (force/area) 

Card  5 Ultimate  stress  (force/area) 

Card  6 Strain  at  which  the  ultimate  stress  is  achieved.  Should 
be  consistent  with  the  equivalent  plastic  strain  defini- 
tion in  Equation  (39)  which  is  essentially  a true 
strain.  Stress  varies  linearly  between  the  yield  stress 
and  the  ultimate  stress.  The  stress  for  larger  strains 
is  equal  to  the  ultimate  until  fracture  occurs. 

Card  7 Maximum  negative  hydrostatic  pressure  (force /area) 

Card  8 Strain  rate  effect  constant,  C^,  in  Equation  (40) 

Card  9 Pressure  effect  constant,  Cg,  in  Equation  (40) 

Card  10  Pressure  effect  constant,  C , in  Equation  (40) 

Card  11  Temperature  effect  constant,  C^,  in  Equation  (40) 
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Card  12  Temperature  effect  constant,  C,.,  in  Equation  (40) 

Card  13  Hydrostatic  pressure  constant,  in  Equation  (42) 
(force/ area) 

Card  14  Hydrostatic  pressure  constant,  in  Equation  (42) 
(force  /area) 

Card  15  Hydrostatic  pressure  constant,  K^,  in  Equation  (42) 
(force/area) 

Card  16  Gruneisen  coefficient,  p,  in  Equation  (42) 

Card  17  Linear  artificial  viscosity  coefficient,  C^  , in  Equa- 
tion (45) 

Card  18  Quadratic  artificial  viscosity  coefficient,  C^,  in 

Q 

Equation  (45) 


Card  19  Equivalent  strain 
Card  20  Volumetric  strain 


If  either  is  exceeded,  there  is  a 
shear  and  tensile  failure.  Posi- 
tive hydrostatic  pressure  and 
viscosity  stress  capability  re- 
main, If  Equivalent  strain  is 
negative,  material  behaves  like 
a liquid. 

Card  21  Equivalent  strain  If  exceeded,  the  element  fails  to- 
tally and  produces  no  stresses  or 
pressures. 

Card  22  Initial  temperature  (degrees) 

Projectile  Scale/Shift/Rotate  Card  (6F10^  0)  - 

XSCALE  = Factor  by  which  the  x coordinates  of  all  projec- 
tile nodes  are  multiplied.  Applied  after  the  co- 
ordinate shifts  (XSHIFT,  ZSHIFT)  described 
later. 
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YSC  ALE 


ZSC  ALE 


XSHIFT 


ZSHIFT 

ROTATE 


Factor  by  which  the  y coordinates  are  multiplied. 

Factor  by  which  the  z coordinates  are  multiplied. 

Increment  added  to  the  x coordinates  of  all  pro- 
jectile nodes  (length).  Applied  before  the  scale 
factors  (XSCALE,  YSC  ALE,  ZSCALE). 

Increment  added  to  the  z coordinates  (length). 

Rotation  about  the  axis  (at  x = z = 0)  of  all  projec- 
tile nodes  (degrees).  Applied  after  the  coordinate 
shifts  (XSHIFT,  ZSHIFT)  and  the  scale  factors 
(XSCALE,  YSCALE,  ZSCALE). 


• Node  Data  Cards  for  Projectile  - These  are  described  in  Section 
IV-l-b.  End  with  a blank  card. 


• Target  Scale/Shift/Rotate  Card  (6F10.  0)  - Same  as  Projectile 
Scale/Shift/Rotate  Card  except  it  applies  to  the  target  nodes. 


• Node  Data  Cards  for  Target  - These  are  described  in  Section 
IV-l-b.  End  with  a blank  card. 


• Element  Data  Cards  for  Projectile  - These  are  described  in 
Section  IV-l-c.  End  with  a blank  card. 


• Element  Data  Cards  for  Target  - These  are  described  in  Section 
IV-l-c.  End  with  a blank  card. 

• Sliding  Surface  Identification  Card  (1115.  5X,  2F10  0)  - Each 
sliding  surface  contains  one  Identification  Card  and  cards  (if 
required)  describing  the  slave  nodes.  For  each  sliding  sur- 
face all  the  slave  nodes  must  have  lower  node  numbers  than 
all  the  master  nodes  of  that  particular  sliding  surface.  The 
slave  nodes  of  a subsequent  sliding  surface  must  all  have 
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higher  node  numbers  than  those  of  the  slave  nodes  of  the  pre- 
vious sliding  surface.  If  there  is  more  than  one  sliding  sur- 
face, all  data  for  the  first  surface  are  input  before  beginning 
data  for  the  second  surface.  The  mass  and  spacing  of  the 
slave  nodes  should  not  be  significantly  greater  than  that  of 
the  master  nodes  in  the  initial  or  deformed  geometry.  Also, 
the  slave  nodes  cannot  be  restrained  in  the  z direction  End 
with  a blank  card  after  the  last  sliding  surface 


TYPE 


1 gives  the  specialized  search  routine  as  de- 
scribed in  Section  II.  The  slave  surface  must 
be  above  the  master  surface  relative  to  the  z 
axis,  and  it  must  be  a single- valued  function  of 
x and  y. 


2 gives  the  generalized  search  routine  as  de- 
scribed in  Section  II. 


IM1 


NML 


NMW 


The  node  number  of  the  first  master  node  which 
corresponds  to  node  1 in  Figure  3. 

Number  of  nodes  per  row  of  master  nodes.  NML 
is  equal  to  N in  Figure  3.  Each  row  of  master 
nodes  must  have  the  same  number  of  nodes. 

Number  of  rows  of  master  nodes.  NMW  is  equal 
to  M in  Figure  3. 

Note:  When  looking  from  the  slave  node  side  to 

the  master  node  side  of  the  sliding  surface, 
each  successive  row  of  master  nodes  is  to 
the  left  of  the  preceding  row  of  nodes  when 
viewed  from  node  1 to  node  N. 
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The  node  number  increment  along  the  rows  of 
master  nodes.  If  IM1  = 100,  NMI,  = 6,  and  IDL  = 
2,  then  the  first  row  of  nodes  in  the  master  sur- 
face consists  of  nodes  100,  102,  104,  10G,  108, 
110. 


1L)W 


IDIA 


IS1 


The  node  number  increment  between  the  first  node 
of  each  row.  If  IDW  = 20  and  IM1,  NML  and  IDL 
are  as  described  in  the  preceding  description  of 
IDL,  then  the  second  row  of  master  nodes  consists 
of  nodes  120,  122,  126,  128,  130. 

1 is  for  the  diagonal  orientation  shown  in  Figure  3. 

/ 2 is  for  the  other  diagonal  orientation  where  the 

diagonals  go  in  the  general  direction  from  the 
first  master  node  to  the  last  master  node. 

The  first  (lowest  node  number)  slave  node. 


ISN  = The  last  (highest  node  number)  slave  node.  Must 

have  a lower  node  number  than  all  the  master  nodes. 

NSLV  = The  total  number  of  slave  nodes  in  the  sliding  sur- 
face. If  all  nodes  between  IS1  and  ISN  are  slave 
nodes  (NSLV  = ISN  - IS1  + 1),  no  additional  input 
data  are  necessary. 


Otherwise,  slave  nodes  are  read  individually  or  in 
groups. 


NSG  = Numbers  of  groups  of  slave  nodes  to  be  read.  If 
NSG  = 0,  the  nodes  are  either  consecutive  between 
IS1  and  ISN,  or  they  are  read  individually. 
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FRICT  = The  coefficient  of  sliding  friction. 

DREF  = A reference  distance  used  to  determine  the  zone 

of  acceptability,  -asshown  in  Figure  3.  Used  only 
for  the  generalized  search  routine  (TYPE  = 2) 

This  distance  should  be  large  enough  such  that  a 
sphere  with  a radius  of  DREF  will  completely 
contain  each  deformed  triangular  master  plane  if 
the  center  of  the  sphere  coincides  with  the  center 
of  the  triangle.  Any  slave  nodes  outside  this 
spherical  zone  will  not  be  considered  for  that  spe- 
cific triangular  plane. 

• Individual  Slave  Node  Cards  (1615)  - Individual  slave  nodes  are 
read  if  the  slave  nodes  are  not  consecutive  between  IS1  and  ISN, 
and  if  the  slave  nodes  are  not  read  in  groups  (NSG  = 0)  Slave 
nodes  must  be  read  in  ascending  order  from  IS1  to  ISN. 

• NSG  Grouped  Slave  Node  Cards  (315)  - 

ISG  1 = The  first  (lowest  number)  node  in  a group  of  nodes 

ISGN  = The  last  (highest  number)  node  in  a group  of  nodes 

INC  = The  increment  between  the  slave  nodes.  If  ISG1  = 

100,  ISGN  = 5,  then  nodes  100,  ’OS,  110,  115,  120 
will  be  designated  as  slave  nodes. 

• Initial  Velocity  Card  (7F10,  0)  - This  card  describes  the  initial 
velocity  conditions  of  the  projectile  and  the  target.  If  there  are 
interface  nodes  which  include  mass  from  both  the  projectile  and 
the  target,  the  velocities  of  these  nodes  are  adjusted  by  the 
program  to  conserve  momenta. 

PXDOT  = Projectile  velocity  in  the  x directum 

(distance/time) 
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Projectile  velocity  in  the  y direction 

Projectile  velocity  in  the  z direction 

Target  velocity  in  the  x direction 

Target  velocity  in  the  y direction 

Target  velocity  in  the  z direction 

Integration  time  increment  (At  in  equations 
of  motion  in  Section  II)  for  the  first  cycle. 

This  must  be  less  than  the  time  required 
to  travel  across  the  minimum  altitude  of 
each  tetrahedron  element  at  the  sound  speed 
of  the  material  in  that  element. 

b.  Input  Data  for  Node  Geometry 

A summary  of  input  data  for  nodal  geometry  is  given  in  Figure  14.  Nodes 
may  be  input  as  a line  of  nodes,  and/or  special  shapes  consisting  of  rods, 
various  node  geometries,  flat  plates  and  spheres.  There  is  no  limit  to 
the  number  of  shapes  included  for  the  projectile  or  the  target.  In  all  cases 
the  nodes  are  numbered  consecutively.  If  a node  is  at  the  interface  of  the 
projectile  and  the  target  and  contains  mass  from  both  the  projectile  and  the 
target,  it  must  be  included  with  the  projectile  nodes. 

Line  of  Nodes  Card  --  Node  cards  for  each  line  of  nodes  are  supplied  as 
needed  by  the  user.  Refer  to  Figure  16  for  the  spacing  of  the  nodes. 

Line  of  Nodes  Description  Card  (215,  6F10.0,  F7.  0,  311)  - 
1 = Identification  number  for  a line  of  nodes. 

NNODE  = Number  of  nodes  in  the  row  of  nodes.  The  nodes 

are  numbered  consecutively. 


PYDOT 

PZDOT 

TXDOT 

TYDOT 

TZDOT 

DTI 
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total  length 
NUMBER  INCREMENTS 


A,  _ N (1  - EXPAND) 

4 


( I - EXfANO*) 


As 

A 


EXPAND' 


' " * EXPAND  1 


Figure  16.  Nodal  Spacing  for  Various  Expansion  Factors 
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XI 


Y1 

Z1 

XN 


YN 

ZN 

EXPAND  = 


IX 


IY 

IZ 


x coordinate  of  first  node  N1  (length) 

y coordinate  of  first  node  N1  (length) 

z coordinate  of  first  node  N1  (length) 

x coordinate  of  last  node  NN  (length).  Leave 
blank  if  a single  node  is  entered. 

y coordinate  of  last  node  NN  (length) 

z coordinate  of  last  node  NN  (length) 

Factor  by  which  the  distance  between  adjacent 
nodes  changes.  For  example,  if  the  distance 
between  the  first  two  nodes  is  A^ , the  distance 
between  the  second  and  third  nodes  is  Ag  = Aj* 
EXPAND.  A summary  of  relative  spacing  be- 
tween the  first  two  nodes  and  the  last  two  nodes 
is  given  in  Figure  16. 

1 will  restrain  all  nodes  (N1  . . NN)  in  the  x di- 
rection. No  restraint  if  left  blank. 

1 will  restrain  nodes  in  y direction. 

1 will  restrain  nodes  in  z direction. 


Rod  Shape  Node  Cards  --  The  rod  geometry  descriptions  are  given  in  Fig- 
ure 17.  The  rod  is  always  generated  in  a vertical  position  about  the  z axix. 
When  viewed  from  the  positive  z direction,  the  nodes  are  numbered  conse- 
cutively counterclockwise,  inner  to  outer  and  downward.  Only  one-half  the 
rod  is  generated  as  shown,  and  restraints  are  provided  normal  to  the  plane 
symmetry.  The  rotation  of  the  rod  for  oblique  impact  is  obtained  with  a 
Scale/Shift/Rotate  Card. 

• Rod  Node  Identification  Card  (15)  - 

2 = Identification  number  for  rod  geometry. 
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NOR 


1 


1 


\ 


Node  Description  Card  (415,  3F10.  0)  - 


NIR 

NPT,N 


IRAD 

ZTOP 

ZBOT 

EXPAND 


Outer  node  ring  number. 

Inner  node  ring  number. 

The  number  of  cross-sectional  planes  of  nodes 
in  the  rod. 

0 gives  ROTOP,  RITOP.  ROBOT,  RIBOT  in- 
put option  for  uniform  radial  spacing. 

1 gives  RT(I),  RB(I ) input  option  (I  = NIR,  NOR). 
The  z coordinate  of  the  top  of  the  rod. 

The  z coordinate  at  the  bottom  of  the  rod. 

Factor  by  which  the  vertical  distance  between 
adjacent  node  changes.  Factor  applies  from 
top  to  bottom.  Same  as  described  for  the  Line 
of  Nodes  Description  Card. 


• Rod  Node  Radii  Card  for  IRAD  = 0 (4F10.  0)  - 


ROTOP 

RITOP 

ROBOT 

RIBOT 


Radius  of  top  outer  node  ring  of  rod. 
Radius  of  top  inner  node  ring  of  rod. 
Radius  of  bottom  outer  node  ring  of  rod. 
Radius  of  bottom  inner  node  ring  of  rod. 


• Rod  Node  Top  Radii  Card  (s)  for  IRAD  = 1 (8F10.0)  - 


RT(I)  = Top  radii  for  each  node  ring,  I = NIR,  NOR. 

Array  RT  is  currently  dimensioned  for  a maxi- 
mum outer  ring  number  of  1 6 (NOR.  LE.  16). 
Mate:  If  NIR  = 0,  then  RT(0)  internally  set  to 
0 and  I = 1,  NOR. 
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Rod  Node  Bottom  Radii  Card  (s)  for  IRAD  - 1 (8F10.  0)  - 


RB(I) 


Bottom  radii  for  each  node  ring,  I = NIR,  NOR. 
Array  RB  is  currently  dimensioned  for  a maxi- 
mum outer  ring  number  of  16  (NOR.  LE.  16). 
Note:  If  NIR  = 0,  then  RB  (0)  internally  set  to 
0 and  I = 1,  NOR. 


Note:  If  it  is  not  possible  to  describe  the  node  geometry  of 
the  rod  with  a single  shape,  it  is  possible  to  use  mul- 
tiple shapes  to  form  a single  rod.  The  nodes  must  be 
numbered  consecutively,  and  the  radii  and  the  number 
of  node  rings  must  be  the  same  for  the  individual  rod 
shapes  at  their  interface.  Also,  ZBOT  and  ZTOP  for 
adjoining  rods  should  not  be  identical.  ZTOP  for  the 
lower  rod  should  be  less  than  ZBOT  for  the  upper  rod 
by  the  desired  node  spacing  in  the  z direction. 


Nose  Shape  Node  Cards  --  The  nose  geometry  descriptions  are  given  in  Fig- 
ure 18.  The  nose  shape  is  always  generated  in  a vertical  position  (pointed 
downward)  about  the  z axis.  When  viewed  from  the  positive  z direction,  the 
nodes  are  numbered  consecutively,  counterclockwise,  downward  and  inner 
to  outer.  Only  one-half  the  nose  shape  is  generated  as  shown,  and  restraints 
are  provided  normal  to  the  plane  of  symmetry.  The  node  geometry  for  the 
flat  plate  at  the  rod  interface  is  not  generated  with  the  nose  generator  and  must 
therefore  be  generated  with  the  rod  generator.  The  number  of  rings  must  be 
identical  for  the  rod  and  the  nose. 


# Nose  Node  Identification  Card  (15)  - 

3 = Identification  number  of  nose  geometry. 

• Nose  Node  Description  Card  (415,  4F10.  0)  - 

NOR  = Outer  node  ring  number. 

NIR  = Inner  node  ring  number. 


i 
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ELEMENT  NODE 

RING  RINC 

5 I 4 | 3 I 2 I 1 I 1 2 3 4 5 


\+—  ROTOP 


ROD 

INTERFACE 


^ ZMIN  

NIR-0  NOR-5  NIER-1  NOER-5  MR-2  NOR-5  NIER-3  NOER-5 

OGIVAL  NOSE  (I  NOSE-3) 


Figure  18.  Geometry  for  Various  Nose  Shapes 


INOSE 


IRAD 


ROTOP 

RITOP 


ZTOP 


ZMIN 


1 identifies  a conical  nose  shape. 

2 identifies  a rounded  nose  shape.  If  the  length 
of  the  nose  is  equal  to  the  radius,  a spherical 
shape  is  generated.  When  the  length  is  not  equal 
to  the  radius,  the  axial  coordinates  are  scaled 
and  the  various  radii  are  not  changed 

3 identified  a tangent  ogival  nose  shape. 

0 gives  ROTOP,  RITOP,  ZTOP,  ZMIN  input 
option  for  uniform  spacing. 

1 gives  RT(I),  ZM(I),  ZTOP  input  option  (I  = 
NIR,  NOR) 

Top  outer  node  radius  of  nose. 

Top  inner  node  radius  of  nose. 

The  z coordinate  at  the  top  of  the  nose.  This  is 
identical  to  ZBOT  for  the  rod  shape  at  the  rod 
interface. 

The  z coordinate  at  the  tip  of  the  nose 


Nose  Node  Top  Radii  Card(s)  for  IRAD  = 1 (8F10.  0)  - 

RT(I)  = Top  radius  for  each  node  ring,  I = NIR,  NOR. 

Array  RT  is  currently  dimensioned  for  a maxi- 
mum outer  ring  number  of  16  (NOR.  LE.  16). 
Note:  If  NIR  = 0,  then  RT(0)  internally  set  to 
0 and  I = 1,  NOR. 

Nose  Node  ZMIN  Card(s)  for  IRAD  = 1 (8F10.  0)  - 


ZM(I) 


Minimum  coordinates  for  each  node  ring,  I = NIR, 
NOR.  Array  ZM  is  currently  dimensioned  for  a 
maximum  outer  ring  number  of  16  (NOR.  LE  16). 
Note:  If  NIR  = 0,  then  ZM(0)  internally  set  to  ZTOP 
and  1 = 1,  NOR. 


Flat-Plate  Shape  Node  Cards  --  The  flat-plate  descriptions  are  given  in 
Figure  19.  In  all  cases,  the  lines  connecting  the  adjacent  corner  nodes  are 
parallel  to  one  of  the  three  primary  axes.  The  nodes  are  generated  in  rows 
parallel  to  the  x axis  and  are  numbered  consecutively  within  each  row  in  the 
direction  of  the  increasing  x axis.  This  generation  continues  row  by  row  in 
the  positive  y direction,  plane  by  plane  in  the  negative  z direction. 

• Flat-Plate  Node  Identification  Card  (15)  - 


4 = Identification  number  for  flat-plate  geometry. 

Flat- Plate  Node  Description  Card  (615,  5F10.  0 


NX  = The  number  of  nodes  in  the  x direction;  in 

Figure  19,  NX  = 13, 

NY  = The  number  of  nodes  in  the  y direction:  in 

Figure  19,  NY  = 7. 

NZ  = The  number  of  nodes  in  the  z direction;  in 

Figure  19,  NZ  = 5. 

NXEND  = The  number  of  nodes  in  the  x direction  in 
each  of  the  two  variable  x spacing  regions. 

The  spacing  is  determined  byX-EXPAND 
and  the  fractional  length  by  X- PART.  In 
Figure  19,  NXEND  = 4.  Depending  on  wheth- 
er NX  is  odd  or  even,  NXEND  can  have  a 
maximum  value  of  either  (NX  +l)/2  or  NX/2, 
respectively,  unless  the  special  option  dis- 
cussed in  X-PART  is  used.  The  remaining 
middle  x region  (if  any)  is  uniformly  spaced. 

NYEND  = The  number  of  nodes  in  the  y direction  in  the 
variable  y spacing  region.  Spacing  is  deter- 
mined by  Y- EXPAND  and  the  fractional  length 


by  Y-PART.  In  Figure  19,  NYEND  = 4. 
NYEND  can  have  a maximum  value  of  NY. 
The  remaining  y region  (if  any)  is  uni- 
formly spaced. 

IY  = 1 gives  restraint  in  the  y direction  only  if 

y = 0 (Y1  * 0.  0). 

X-E'XPAND  - Factor  by  which  the  x distance  between  ad- 
jacent nodes  changes  outward  to  the  ends  for 
each  X-PART  variable  spacing  region. 

Same  as  described  for  the  Line  of  Nodes 
Description  Card 

X-PART  = Fractional  part  of  the  total  x length  of  the 
flat  plate  occupied  by  each  of  the  varia- 
ble x spacing  regions. 

Note:  If  X-PART  = 0.  0,  the  entire  spacing 
in  the  x direction  is  uniform.  If  X-PART  = 

1.  0 and  NXEND  = NX,  the  spacing  on  the 
positive  x direction  is  variable  for  entire  x 
length  (XI  to  XN) 

Y-EXPAND  = Factor  by  which  the  y distance  between  ad- 
jacent nodes  changes  in  the  increasing  y 
direction  for  the  Y-PART  variable  spacing 
region. 

Y - PART  = Fractional  part  of  the  total  y length  of  the 

flat  plate  occupied  by  the  variable  y spacing 
region. 

Z-EXPAND  = Factor  by  which  the  z distance  between  ad- 
jacent nodes  changes  in  the  decreasing  z 
direction  from  Z1  to  ZN. 
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• Flate-Plate  Node  Size  Card  (6  F 10.0)  - 

XI  = The  minimum  x coordinate  of  the  plate  shape 

Y1  = The  minimum  y coordinate  of  the  plate  shape 

Z1  = The  maximum  z coordinate  of  the  plate  shape 

XN  = The  maximum  x coordinate  of  the  plate  shape 

YN  = The  maximum  y coordinate  of  the  plate  shape 

ZN  = The  minimum  z coordinate  of  the  plate  shape 

Sphere  Shape  Node  Cards  --  The  cross-section  of  the  bottom  spherical 
shape  is  identical  to  that  shown  for  the  rounded  nose  shape  of  Figure  18. 

The  top  one-half  cross-section  is  initially  geometrically  summetric  to  the 
bottom.  Only  one-half  of  the  top  and  bottom  halves  are  generated,  and  re- 
straints are  provided  normal  to  the  vertical  plane  of  symmetry.  The  sphere 
is  generated  with  the  nodes  numbered  as  two  rounded  circular  noses  having 
an  interface  between.  The  top  nose  is  generated  first;  viewed  from  the 
positive  z direction,  this  generation  is  counterclockwise,  upwards  and  inner 
to  outer.  The  bottom  nose  is  generated  with  the  interface  included  with  each 
spherical  shell;  this  generation  viewed  from  the  positive  z direction  is  counter- 
clockwise, downwards  and  inner  to  outer.  A summary  of  the  number  of  nodes 
included  in  the  various  shapes  is  given  in  Figure  20. 

• Sphere  Node  Identification  Card  (15)  - 


5 = Identification  number  for  sphere  geometry. 

Sphere  Node  Description  Card  (215,  5X,  115,  3F10  0)- 

NOR 

= 

Outer  node  ring  number. 

NIR 

= 

Inner  node  ring  number. 

IRAD 

f 0 gives  RO,  RI,  ZCG  input  option  for  uniform  spacing 

1 1 gives  RT(I),  ZCG  input  option  (I  = NIR. 

^ NOR). 
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DOES  NOT  INCLUDE  NODES  AT  ROD  INTERFACE 


Figure  20.  Summary  of  Nodes  and  Elements  in  Various  Shapes 


RO  = Radius  of  outer  node  ring. 

RI  = Radius  of  inner  node  ring. 

ZCG  = The  z coordinate  of  the  center  of  the  sphere. 

• Sphere  Node  Radii  Card(s)  for  IRAD  = 1 (8F10,  0)  - 

RT(I)  = Radii  for  each  node  ring,  I = KIR,  NOR. 

Array  RT  is  currently  dimensioned  for  a 
maximum  outer  ring  number  of  16  (NOR. 

LE.  16). 

Note:  If  NIR  = 0,  then  RT(0)  internally  set 
to  0 and  1 = 1,  NOR. 


c.  Input  Data  for  Element  Geometry 

A summary  of  input  data  for  element  geometry  is  given  in  Figure  15.  Ele- 
ments maybe  input  as  a series  of  individual  or  composite  elements  and/or 
special  shapes  of  rods,  nose  geometries,  flat  plates,  and/or  spheres.  The 
elements  must  be  assembled  in  a manner  consistent  with  the  previously  gen- 
erated nodal  geometry.  There  is  no  limit  to  the  number  of  shapes  included 
for  the  projectile  or  the  target. 

Series  of  Composite  Elements  Card  --  Element  cards  for  each  Series  of 
Composite  Elements  are  supplied  as  needed  by  the  user.  For  this 
discussion  it  will  be  assumed  that  the  elements  are  entered  as  a series  of 
composite  brick  elements,  each  containing  six  individual  elements  as  shown 
in  Figure  21.  Following  this  immediate  discussion  will  be  an  example  and 
instructions  for  generating  a series  of  individual  elements. 

• Series  of  Composite  Elements  Description  Card  (1215)  - 

1 = Identification  number  for  a series  of  composite 

elements. 
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NCOMP 

MATL 


N1-N8 


INC 


Number  of  composite  elements  in  the  series. 

Material  number  (1,  2,  3,  4,  or  5)  of  the 
elements.  If  left  blank,  the  material  number 
from  the  previous  element  data  card  will  be 
used. 

Node  numbers  of  the  first  composite  brick 
element  as  shown  in  Figure  2b,  Nodes  Nl, 
N2,  N3,  N4  and  nodes  N5,  N6,  N7,  N8  are 
counterclockwise  when  looking  from  Nl  to 
N 5 . 

The  node  number  increment  added  to  the  node 
numbers  of  the  previous  composite  brick  ele- 
ment for  the  next  composite  brick  element. 


An  example  of  input  data  for  composite  brick  elements  is  shown  in  Figure  22. 
In  the  upper  left  it  can  be  seen  that  there  are  four  rows  of  nodes  (1  to  4,  5 to 

8,  9 to  12,  13  to  16)  which  are  arranged  to  contain  three  composite  brick 
elements.  If  the  first  element  is  numbered  1,  then  the  first  composite  brick 
contains  elements  1 to  6,  the  second  contains  7 to  12,  and  the  third  contains 
13  to  18.  The  first  composite  brick  is  defined  by  nodes  Nl  =1,  N2  = 5,  N3  = 

9,  N4  = 13,  N5  =2,  N6  = 6,  N7  = 10,  and  N8  = 14.  Note  that  Nl  to  N4  and  N5 
to  N8  are  counter-clockwise  when  looking  from  Nl  to  N5.  The  six  individual 
elements  are  generated  according  to  the  arrangement  and  order  (A,  B,  C,  D, 
E,  F)  shown  in  Figure  21.  The  node  numbers  for  each  successive  brick  are 
simply  INC  = 1 greater  than  those  of  previous  brick.  For  the  second  brick  for 
instance,  Nl  = 1 + 1 = 2,  N2  = 5 + 1 = 6,  N3  = 9 + 1 = 10,  N4  = 13  + 1 = 14, 

N5  = 2 + 1 = 3,  N6  = 6 + 1 = 7,  N7  = 10  + 1 = 11,  and  N8  = 14  + 1 = 15. 


It  is  possible  to  generate  a series  of  individual  tetrahedron  elements  by  letting 
Nl  to  N4  be  the  nodes  of  the  first  element,  where  Nl,  N2,  and  N3  are  counter- 
clockwise when  viewed  from  N4.  This  option  is  exercised  when  N5  to  N8  are 
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N2- 5 N3-  9 

BRICK  1 
ELEMENTS  I -6 


BRICK  2 

ELEMENTS  7-12 


Figure  22.  Node  Element  Input  Data  Exampl 


left  blank.  It  is  also  possible  to  generate  a series  of  composite  wedge  ele- 
ments, each  containing  three  individual  tetrahedron  elements.  The  three 
elements  in  a composite  wedge  element  are  numbered  consecutively.  If  N2 
and  N6  are  left  blank,  the  first  three  tetrahedron  elements  (A,  B,  C)  are  de- 
fined by  nodes  Nl,  N3,  N4,  N5,  N7  and  N8  as  shown  in  Figure  21,  Likewise, 
if  N4  and  N8  are  left  blank,  the  first  three  elements  (D,  E,  F)  are  defined  by 
nodes  Nl,  N2,  N3,  N5,  N6,  and  N7. 

Rod  Shape  Element  Cards  --  Elements  are  generated  for  the  rod  shapes  illus- 
trated in  Figure  17  and  described  by  the  Rod  Shape  Node  Cards.  The  elements 
are  numbered  consecutively  and  are  generated  in  layers  of  composite  brick 
elements  beginning  with  top  layer  1 and  ending  at  bottom  layer  NLA  Y.  The 
entire  first  layer  of  elements  is  generated  before  the  second  layer,  etc.,  and 
the  composite  brick  elements  of  each  layer  are  generated  in  a counterclock- 
wise manner  for  each  ring  of  elements  from  the  inner  to  the  outer  ring. 

• Rod  Element  Identification  Card  (15)  - 

2 = Identification  number  for  rod  geometry, 

• Rod  Element  Description  Card  (615)- 

NOER  = Outer  element  ring  number, 

NIER  = Inner  element  ring  number, 

NLAY  = The  number  of  layers  of  elements  in  the  rod. 

The  total  number  of  elements  in  a rod  shape 
shown  in  Figure  17  is  dependent  on  the  num- 
ber of  layers  and  the  number  elements  per 
layer.  The  number  elements  per  layer  is  de- 
pendent on  the  inner  and  outer  element  ring 
numbers.  [For  example:  If  NLAY  = 10, 

NIER  = 2 and  NOER  = 5,  the  total  number  of 
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IMAT 


MATL 


elements  can  be  determined  through  use  of 
Figure  2 0.  The  number  of  elements  per 
layer  for  the  solid  rod  of  NOFIR  = 5 is  300, 
and  of  NIER  = 2 is  48  Therefore  the  total 
of  elements  for  the  hollow  cylinder  is  10* 
(300-48)  = 2520] 

The  number  of  the  lowest  numbered  rod  node. 
For  the  solid  rod,  this  is  the  centerline  node 
on  the  top  end  of  the  rod.  For  the  hollow  rod, 
this  is  the  innermost  clockwise  node  on  the 
top  end  of  the  rod  when  viewed  from  the  top. 

0 gives  MATL  input  option.  Material  num- 
ber of  all  element  rings  uniformly  set  to 
MATL. 

1 gives  M(I)  input  option  (I  = NIER,  NOER). 
Material  numbers  for  each  element  ring  must 
be  input. 

Material  number  (1.  2,  3,  4 or  5)  of  the  uni- 
form material  rod. 


• Rod  Element  Card  for  IMAT  = 1 (1615)  - 


M(I)  = Material  number  for  each  element  ring,  I = 

NIER,  NOER.  Array  M is  currently  dimen- 
sioned for  a maximum  outer  ring  number  of 
16  (NOER.  LE,  16), 

Nose  Element  Cards  --  Elements  are  generated  for  the  nose  shapes  illustrated 
in  Figure  18  and  described  by  the  Nose  Shape  Node  Cards.  The  elements  are 
numbered  consecutively  and  are  generated  in  shells  of  composite  brick  ele- 
ments beginning  with  the  innermost  shell  and  ending  with  the  outermost  shell. 
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The  entire  first  shell  of  elements  is  generated  before  the  second  shell,  etc.  , 
and  the  composite  brick  elements  of  each  shell  are  generated  in  a counter- 
clockwise manner  for  each  ring  of  elements  from  the  top  to  the  bottom  of 
each  shell. 


• Nose  Element  Identification  C ard  (15)- 

3 = Identification  number  for  nose  geometry. 


Nose  Element  Description  Card  (215.  5X.  315) 


NOER 


Outer  element  ring  number. 


NIER 


Inner  element  ring  number. 


N1  = The  number  of  the  lowest  numbered  nose 

node.  (The  nose  does  not  include  interface 
nodes, ) 


IMAT 


MATE 


0 gives  MATL  input  option.  Material  num- 
ber of  all  element  rings  uniformly  set  to 
MATL. 


1 gives  M(I)  input  option  (I  = NIER,  NOER). 
Material  number  for  each  element  ring  must 
be  input. 


Material  number  (1,  2,  3,  4,  or  5)  of  the  uni- 
form material  nose. 


>se  Element  Material  Card  for  IMAT  = 1 


M(I)  = Material  number  for  each  element  ring,  I = 

NIER,  NOER.  Array  M is  currently  dimen- 
sioned for  a maximum  outer  ring  number  of 
16  (NOER.  LE.  16). 
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Flat- Plate  Shape  Element  Cards  --  The  elements  are  generated  for  the  flat- 
plate  illustrated  in  Figure  19  and  described  by  the  FI  at- Plate  Shape  Node 
Cards.  The  elements  are  numbered  consecutively  and  are  generated  in  rows 
of  composite  brick  elements.  The  rows  of  elements  go  in  the  direction  of 
the  increasing  x axis.  Each  row  of  elements  is  to  the  positive  y direction  of 
the  preceding  row.  After  a layer  is  complete,  the  next  lower  layer  is  gener- 
ated in  a similar  manner. 

• Flat-Plate  Element  Identification  Card  ( 15 ) ~ 

4 = Identification  number  for  flat- plate  geometry 

• Flat- Plate  Element  Description  Card  (415.  5X.  15)- 

NLX  = Number  of  layers  of  composite  brick  elements 

in  the  x direction  The  total  number  of  nodes 
along  the  x direction  must  be  NLX  + 1.  In 
Figure  19,  NLX  = 12. 


NLY 

NLZ 

N1 

MATL 


Number  of  layers  of  composite  brick  elements 
in  the  y direction.  The  total  number  of  nodes 
along  the  y direction  must  be  NLY  +1.  In 
Figure  19,  NLY  * 6. 

Number  of  layers  of  composite  brick  elements 
in  the  z direction.  In  Figure  19,  NLZ  = 4 

The  node  number  of  the  corner  node  shown 
in  Figure  19. 

Material  number  (1,  2,  3,  4 or  5)  of  the  flat- 
plate. 


Sphere  Shape  Element  Cards  --  Elements  are  generated  for  a sphere,  the  bot- 
tom half  cross  section  of  which  is  identical  to  the  rounded  nose  shown  in  Fig- 
ure 18  and  described  by  the  Sphere  Shape  Nose  Cards  When  viewed  from  the 
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top,  the  elements  are  consecutively  numbered  counterclockwise,  upwards  and 
outwards  for  the  top  one-half  and  then  counterclockwise,  downwards  and  out- 
wards for  the  bottom  one-half. 


• Sphere  Element  Identification  Card  (15)- 

5 = Identification  number  for  sphere  shape. 


• Sphere  Element  Description  Card  (215,  5X,  315)- 


NOER 
NIER 
N 1 


Outer  element  ring  number. 

Inner  element  ring  number. 

The  number  of  the  lowest  numbered  sphere 
node 


IMAT 


MATL 


r 0 gives  MATL  input  option.  Material  num- 
ber of  all  element  rings  uniformly  set  to 
MATL. 

1 gives  M(I)  input  option  (I  = NIER,  NOER). 
Material  number  for  each  element  ring  must 
v be  input. 


= Material  number  (1,  2,  3,  4 or  5)  of  the  uni- 
form material  sphere. 


• Sphere  Element  Material  Card  for  IMAT  * 1 (1615)- 

M(I)  = Material  number  for  each  element  ring,  1 = 

NIER,  NOER.  A rray  M is  currently  dimen- 
sioned for  a maximum  outer  ring  number  of 
16  (NOER.  LE.  16). 
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d.  Sample  Input  Data  for  the  Preprocessor 

Input  data  for  a typical  problem  involving  oblique  impact  of  a copper  rod 
onto  a steel  plate  are  given  in  Section  V. 


2.  INPUT  DATA  FOR  THE  MAIN  ROUTINE 

The  function  of  the  Main  Routine  is  to  perform  the  computations.  The  Main 
Routine  reads  initial  conditions  from  the  restart  tape  which  must  be  generated 
from  a Preprocessor  run  or  a previous  Main  Routine  run.  The  descriptions 
which  follow  are  for  the  data  in  Figure  23.  Consistent  units  must  be  used. 

• Description  Card  - A description  of  the  problem  provided  by 
the  user. 

• Identification  Card  (315.  5X.  F10.  0)  - 
CASE 
CYCLE 


IPRES 


= Case  number  for  run  identification 

= The  cycle  number  at  which  the  restart  occurs. 
The  cycle  numbers  for  which  restart  files  are 
written  are  given  in  the  printed  output  of  the 
previous  run.  (Preprocessor  or  Main  Routine) 

s 0 gives  no  applied  pressures  read  or  applied 

1 will  use  the  pressure  data  which  was  input 
= c in  a previous  run. 

2 will  read  applied  pressures  to  be  used  in 
v subsequent  computations. 
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CPMAX 


! 


EMAX 


Central  Processor  time  (hours)  at  which  the 
results  will  be  written  onto  the  restart  tape 
and  the  run  will  stop.  This  feature  will  be 
bypassed  if  CPMAX  = 0. 

The  upper  limit  for  total  kinetic  energy  if 
applied  pressures  are  included  (IPRES  --  1 or 
2).  This  is  used  for  numerical  instability 
checks.  Run  will  stop  if  the  kinetic  energy 
exceeds  EMAX.  Reave  blank  if  there  are  no 
applied  pressures. 


• Integration  Time  Increment  Card  (4F10.  0)- 


DTMAX  = 

The  maximum  integration  time  increment 
which  will  be  used  for  the  equations  of  mo- 
tion. If  At  from  Equation  (54)  is  greater  than 
DTMAX,  it  will  be  redefined  as  At  = DTMAX. 

DTMIN 

The  minimum  integration  time  increment  al- 
lowed. If  At  from  Equation  (54)  is  less  than 

DTMIN,  the  results  will  be  written  onto  the 
restart  tape,  and  the  run  will  stop. 

SSF 

The  fraction  of  the  sound  speed  transit  time 
used  for  the  integration  time  increment.  This 
is  identical  to  Ct  in  Equation  (54)  and  must  be 
less  than  unity. 

TMAX 

The  maximum  time  the  problem  is  allowed  to 
run.  This  time  refers  to  the  dynamic  response 
of  the  system,  not  the  central  processor  time 
(CPMAX)  described  in  the  Identification  Card. 

The  results  at  time  = TMAX  are  written  onto 
the  restart  tape,  ana  the  run  is  discontinued. 
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Pressure  Cards  for  IPRES  = 2 (615,  FlO.O)-  These  cards 
describe  the  applied  pressures  and  the  elements  to  which 
they  are  applied.  If  other  pressures  were  used  previously 
they  are  all  deleted,  and  the  only  applied  pressures  which 
act  are  those  which  are  input  in  the  current  run.  End  with 
a blank  card. 

ELE1  = The  first  element,  in  a series  of  elements, 
to  which  pressure  is  applied.  Must  not  be 
less  than  ELE1  or  ELEN  from  a previous 
Pressure  Card. 

ELEN  ■ The  last  element  in  series  of  elements. 
Cannot  be  less  than  ELE1. 


ELE  INC 


N1 


NN 


NODE 

INC 


PRES 


The  element  number  increment  between 
ELE1  and  ELEN.  If  ELE1  = 100,  ELEN  = 

12  0 and  ELE  INC  = 5,  then  pressures  are 
applied  to  elements  100,  105,  110,  120. 

The  node  number  opposite  the  triangular 
face  of  element  ELE1  to  which  the  pressure 
is  applied. 

The  node  number  opposite  the  triangular 
face  of  element  ELEN  to  which  the  pressure 
is  applied. 

The  node  number  increment  between  N1  and 
NN.  For  the  elements  described  under  ELE 
INC  (100,  105,  110,  115.  120),  if  N1  = 200, 
NN  = 208  and  NODE  INC  = 2,  then  the  pres- 
sures are  applied  to  the  triangular  faces  op- 
posite nodes,  200,  202,  204,  206,  208,  of 
elements  100,  105,  110,  115,  120. 

The  pressures  which  are  applied  to  the  trian- 
gular faces  of  the  elements  described  on  this 
card  (force/area). 


i 
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• Time-Pressure  Cards  for  IPRES  = 2 (2F10.  0)  - These  cards 
allow  the  applied  pressures  to  be  varied  as  a function  of  time 
A minimum  of  two  cards  must  be  used  which  span  the  time 
from  the  beginning  of  the  run  to  TMAX.  Program  currently 
dimensioned  for  a maximum  of  50  cards.  End  with  a blank 
card. 

PTIME  = The  time  corresponding  to  P(T).  Cards 
must  be  input  in  order  of  increasing 
time . 

P(T)  = The  factor  by  which  all  pressures  are  multi- 

plied at  the  corresponding  time.  Equal  to 
Ft  in  Equations  (48)  to  (50).  Intermediate 
values  are  linearly  interpolated  between 
values  at  specified  times. 

• Data  Output  Cards  (2F10.  0.  15)  - These  cards  are  used  to 
specify  various  forms  of  output  data  at  selected  times,  and 
the  last  card  must  be  for  a time  greater  than  TMAX,  even 
though  output  will  not  be  provided  for  that  specific  time.  Re- 
call that  output  is  automatically  provided  at  TMAX  and  a data 
output  card  need  not  be  provided  for  this  time.  End  run  with 
two  blank  cards. 

TIME  ='  Time  at  which  output  will  be  provided. 

ECHECK  - A code  which  governs  the  printed  output.  Two 
options  are  provided 

(1)  If  ECHECK  is  greater  than  1000,  the  in- 
dividual node  and  element  data  will  not  be 
printed.  Only  system  data  such  as  centers 
of  gravity,  energies,  momenta,  and  net 
velocities  are  provided  for  the  projectile, 
target,  and  total  system. 
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(2)  If  ECHECK  is  less  than  1000,  the  sys- 
tem data  and  individual  node  data  will 
be  printed.  Individual  element  data  will 
be  printed  for  all  elements  which  have 
an  equivalent  plastic  strain  [Equation 
(39)]  equal  to  or  greater  than  ECHECK. 
For  example,  if  ECHECK  = 0,  all  ele- 
ment data  will  be  printed  If  ECHECK  = 
0.  5,  only  those  elements  with  equiva- 
lent plastic  strains  equal  to  or  greater 
than  0.  5 will  have  data  printed.  If 
ECHECK  = 999,  no  element  data  will 
be  printed. 

Note:  There  are  some  instances  when  data 
are  printed  even  though  not  specified 
with  a data  output  card.  If  the  minimum 
time  increment,  DTMIN,  is  violated,  or 
the  specified  central  processor  time, 
CPMAX,  is  exceeded  or  the  maximum 
run  time,  TMAX,  is  achieved,  the  re- 
sults are  written  on  the  restart  tape. 

The  value  of  ECHECK  used  for  the 
printed  output  is  that  value  which  is  on 
the  following  data  output  card 

f 0 will  not  write  results  on  the  restart  tape. 

^ 1 will  write  results  on  the,  restart  tape. 


3.  INPUT  DATA  FOR  THE  POSTPROCESSOR 

The  function  of  the  Postprocessor  is  to  provide  plots  of  the  results  in  the 
form  of  geometry  and  velocity  vectors.  The  Postprocessor  reads  data 
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from  the  restart  tape  which  must  be  generated  by  a Preprocessor  run  or  a 
Main  Routine  run.  The  descriptions  which  follow  are  for  the  data  in  Figure 
24.  Each  plot  requires  a Plot  Identification  Card  and  a Plot  Limits  Card. 
The  three-dimensional  plots  also  require  a 3D  Perspective  Card.  The  data 
must  be  input  in  groups  of  two  (2D  plots)  or  three  (3D  plots)  cards  and  the 
requested  cycle  numbers  must  be  equal  or  greater  than  the  previously  re- 
quested cycle  number.  Consistent  units  must  be  used.  End  with  a blank 
card. 

• Plot  Identification  Card  (315,  5X,  F10.  0,  5A6)  - 


TYPE 

VIEW 

CYCLE 

VSCALE 

TITLE 


1 gives  a geometry  plot 

2 gives  a velocity  vector  plot 

A code  to  specify  the  view  requested.  Op- 
tions are  given  in  Figure  2 6. 

The  cycle  number  of  the  plot  which  is  de- 
sired. The  cycle  numbers  of  the  data 
written  on  the  restart  tape  are  given  in  the 
printed  output  of  the  Preprocessor  and  the 
Ma  in  Routine. 

The  velocity  which  will  give  a velocity  vec- 
tor which  has  a length  of  1.  0 inch  using  the 
scale  of  the  plot.  Leave  blank  for  Geometry 
Plot. 

A title  which  is  written  on  the  plot 


• Plot  Limits  Card  (6F10.  0)  - This  card  specifies  the  portion  of 
the  problem  which  is  plotted  Regions  beyond  those  specified 
are  not  plotted.  For  the  two-dimensional  plots,  the  vertical 
axis  is  10  inches  long  and  the  horizontal  axis  is  as  specified. 


L 


no 


PLOT  IDENTIFICATION  CARD  (31!),  SX  F’U  U.  5A0» 


| TYPE  | VIEVW  |CYCLE|  | VSCAtf  | 


PLOT  LIMITS  CARD  '6F!U  0> 


I 


zmax 


3D  PERSPECTIVE  CARD  FOK  VIEW  • 4 (6F'0  U IS) 

| XPlAHt  | YPLATie  | ZPLANE  |lHID£ 


T 


blank  card 

^ ENDS  RUN 

VIEW 

VIEW 

• Z 

X Z PLOT 

Y Z PLOT 

type  • 1 

1 GEOMETRY  PLOT 

VIEW  ' 

• 3 

X Y PLOT 

TYPE  • 2 

? VELOCITY  PLOT 

VIEW 

■ 4 

3D  PERSPECTIVE  PLOT 

Figure  24.  Postprocessor  Input  Data 


The  length  of  the  horizontal  axis  will  vary  since  it  has  the 
same  scale  as  the  vertical  axis.  The  vertical  axes  are  the 
z axis  (VIEW  = 1 and  2)  and  the  y axis  (VIEW  = 3).  For  the 
three-dimensional  plots,  the  axes  are  scaled  such  that  the 
entire  region  is  included  within  the  10- inch  vertical  axis. 

XMAX 

The  maximum  x coordinate  included  in  the 
plot  (length) 

XMIN 

The  minimum  x coordinate  included  in  the 
plot 

YMAX 

The  maximum  y coordinate  included  in  the 
plot 

YMIN 

The  minimum  y coordinate  included  in  the 
plot 

ZMAX 

The  maximum  z coordinate  included  in  the 
plot 

ZMIN 

The  minimum  z coordinate  included  in  the 
plot 

• 3D  Perspective  Card  for  VIEW  = 4 (6F10.  0,  15)  - 
This  card  is  included  only  for  the  3D  plots. 


XEYE 

YEYE 

ZEYE 

XPLANE 
YPLANE 
Z PLANE 


= ^ Coordinates  of  the  observer  (length) 


Coordinates  included  in  the  plane  on  which 

jthe  results  are  plotted.  The  plane  is  nor- 
mal to  a line  from  XEYE,  YEYE,  ZEYE  to 
XPLANE,  YPLANE,  ZPLANE. 


0 will  plot  alt  free  surfaces  (no  hidden  lines) 

1 will  plot  only  surfaces  which  are  visible  to 
the  observer 
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4.  OUTPUT  DATA 


The  output  data  are  generally  described  using  the  terminology  of  this  report. 

A summary  of  the  printed  output  for  a Preprocessor  run  follows: 

• All  input  data  are  printed  exactly  as  read. 

• Material  data  are  printed  for  the  input  data. 

• All  the  geometric  input  data  are  printed  in  a form  similar 
to  that  used  to  read  the  data. 

• For  each  node,  the  initial  coordinates,  mass 

are  printed  if  IPRINT  = 1 on  the  Identification  Card.  The 
slave  nodes  are  identified  by  adding  "1"  before  the  XYZ  re- 
straint. If  the  XYZ  restraint  is  100,  the  node  is  restrained 
in  the  x direction.  If  the  XYZ  restraint  is  1100,  it  is  a slave 
node  with  the  same  restraint. 

• For  each  element,  the  four  nodes,  volume,  and  material  num- 
ber are  printed  if  IPRINT  = 1. 

• A summary  of  the  geometry,  bandwidth  requirements,  mass 
distributions,  and  initial  dynamic  conditions  is  provided. 

Note:  The  node  and  element  block  sizes  must  always  be  equal 
to  those  used  in  the  Preprocessor.  The  number  of  node 
blocks  can  be  varied,  however,  provided  there  are  at  least 
as  many  as  specified  under  the  bandwidth  requirements. 

• A summary  of  system  data  is  provided  to  give  centers  of 
gravity,  energies,  momenta,  and  net  velocities  for  the  ini- 
tial condition. 

For  a Main  Routine  run,  the  printed  output  data  conforms  to  the  following: 

• All  input  data  are  printed  exactly  as  read. 
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• Pressure  data  are  printed  in  a form  similar  to  that  used  to 
read  the  data 

• For  each  integration  time  increment,  selected  data  are  printed. 
These  consist  of  cycle  number  (CYCLE),  the  time  (TIME),  the 
integration  time  increment  for  the  next  cycle  (DTNEXT),  the 
element  number  which  governs  the  time  increment  (ECRIT), 
and  the  total  kinetic  energy  in  the  system  (TOTAL  K.  E.  ). 

If  ECRIT  = 0 the  integration  time  increment  determined  from 
Equation  (54)  is  either  greater  than  DTMAX,  or  greater  than 
1.  1 times  the  integration  time  increment  for  the  previous 
cycle. 

• Node  data  are  printed  at  the  selected  times  specified  on  the 
Data  Output  Cards  if  ECHECK  is  less  than  1000.  These  data 
include  the  node  number  (NODE),  the  coordinates  (X,  Y,  Z), 
the  velocities  (XDOT,  YDOT,  ZDOT)  and  the  accelerations 
(XDD,  YDD,  ZDD).  If  a slave  node  is  in  contact  with  the  mas- 
ter surface,  or  if  a master  node  is  affected  by  a slave  node, 
then  the  printed  accelerations  are  not  correct  for  these  nodes. 

• Element  data  may  be  printed  at  the  selected  times  specified 
on  the  Data  Output  Cards  if  ECHECK  is  less  than  1000.  Only 
those  elements  with  an  equivalent  plastic  strain  equal  to  or 
greater  than  ECHECK  will  have  their  data  printed.  These  data 
include  the  element  number  (ELE),  the  volumetric  strain 
(DVOL),  the  equivalent  plastic  strain  (EPBAR)  of  Equation 
(39),  the  Von  Mises  equivalent  stress  (VMISES)  of  Equation 
(31),  the  normal  deviator  stresses  (SX,  SY,  SZ)  of  Equations 

( 1 9)  to  (21),  the  shear  stresses  (SXY,  SXZ,  SYZ)  of  Equa- 
tions (22)to(24),  the  net  pressure  (P  + Q)  containing  both  the 
hydrostatic  pressure  of  Equation  (42)  and  the  artificial  vis- 
cosity of  Equation  (45),  and  the  temperature  (TEMP)  of 
Equation  (41). 
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• A summary  of  system  data  is  printed  at  the  selected  times 
specified  on  the  Data  Output  Cards.  Included  are  centers 
of  gravity,  energies,  momenta,  and  net  velocities. 

For  the  Postprocessor,  all  input  data  are  printed  exactly  as  read. 

5.  INSTRUCTIONS  FOR  CHANGING  PROGRAM  DIMENSIONS 

It  may  sometimes  be  desirable  to  change  the  dimensions  of  the  program. 

The  user  should  review  Section  III  to  understand  the  data  structure  of  the 

program. 

To  change  the  size  of  the  node  blocks,  the  following  changes  are  required: 

• Set  NBSIZE  to  the  node  block  size  in  the  data  statements  in 
the  Preprocessor  (PREP),  Main  Routine  (MAIN),  and  Post- 
processor (POST).  This  must  be  consistent  with  the  dimen- 
sion of  the  COMMON/NODE/arravs  (X.  . .MASTER),  and  the 
number  of  node  blocks  in  core,  such  that 

NBSIZE  = NASIZE/NBAND 

where  NBAND  is  the  number  of  node  blocks  in  core  and 
NASIZE  is  the  dimension  of  the  COMMON/NODE/arrays. 

NBAND  and  NASIZE  can  also  be  altered  in  the  same  data 
statements. 

• Dimension  INBUFW  to  (NBSIZE  + 2)  and  RNBUFW  to 
(llaNBSIZE)  in  COMMON/NBUFW/. 

• Dimension  INBUFR  to  (NBSIZE  + 2)  and  RNBUFR  to 
(llaNBSIZE)  in  COMMON/NBUFR/ . 

• Note:  The  node  block  size  cannot  be  changed  during  the 
course  of  a run. 


95 


To  change  the  size  of  the  node  arrays,  the  following  changes  are  required: 

• Dimension  the  COMMON/NODE/arrays  (X. . . MASTER) 
to  the  desired  size. 

• Set  NASIZE  to  the  dimensions  of  the  node  arrays  in  the  data 
statements  in  the  Preprocessor  (PREP),  Main  Routine 
(MAIN),  and  Postprocessor  (POST).  Recall  that 

NASIZE  = NBSIZE  • NBAND 
for  consistency  with  the  node  block  sizes. 

• Note:  The  node  array  size  (NASIZE)  can  be  changed  during 
the  course  of  a run  provided  the  node  block  size  (NBSIZE) 
does  not  change. 

• Note:  Program  is  currently  limited  to  100  node  blocks 
(NBAND  = 100).  It  can  be  increased  by  changing  dimensions 
of  NBTBLE  and  ISLIDE  in  COMMON /BOOKKP/. 

To  change  the  size  of  the  element  block,  the  following  changes  are  required 

• Set  LBSIZE  to  the  element  block  size  in  the  data  statements 
in  the  Preprocessor  (PREP),  Main  Routine  (MAIN)  and  Post- 
processor (POST). 

• Dimension  the  COMMON/ ELEMENT/arra.ys  (NODE1 . . . 

PRES4)  to  the  element  block  size  (LBSIZE) 

• Dimension  IEBUFW  to  (6aLBSIZE  + 2)  and  REBUFW  to 
( 1 5#LBSIZE)  in  COMMON/EBUFW/. 

• Dimension  IEBUFR  to  (6*LBSIZE  + 2)  and  REBUFW  to 
( 1 5*LBSIZE)  in  COMMON/EBUFR/ . 


96 


• Note:  The  element  block  size  (LBSIZE)  cannot  be  changed 
during  the  course  of  a run 

6.  FILE  DESIGNATIONS 

The  files  are  defined  in  the  data  statements  in  the  Preprocessor  (PREP), 
Main  Routine  (MAIN)  and  Postprocessor  (POST).  The  current  designations 
and  descriptions  are: 


NFILER 

= 1 

One  of  two  internal  node  files 

NFILEW 

= 2 

One  of  two  internal  node  files 

LFII.ER 

= 3 

One  of  two  internal  element  files 

LFILEW 

= 4 

One  of  two  internal  element  files 

IN 

= 5 

Input  file 

IOUT 

= 6 

Output  file 

ISFILE 

= 8 

Internal  sliding  surface  file 

ITAPIN 

= 9 

Restart  tape  read  by  MAIN  and  POST 

ITAFOT 

= 10 

Restart  tape  generated  by  PREP  and  MAIN 

ITRIAN 

= 11 

Triangle  data  record  generated  by  POST 

I TEMP 

= 12 

Scratch  file  used  by  POST 

7.  CENTRAL  PROCESSOR  TIMF  ESTIMATES 

For  various  problems  run  on  a CDC  6600  computer,  the  EPIC-3  Main  Rou- 
tine has  used  from  0.  006  to  0.  009  central  processor  second  per  cycle 
per  node. 
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SECTION  V 
EXAMPLES 


Several  examples  have  been  run  to  verify  and  demonstrate  the  capability  of 
the  EPIC -3  code.  Two  wave  propagation  solutions  are  shown  in  Figure  25. 

The  water  column  example  is  for  impact  onto  a rigid  boundary  at  300  meters/ 
second  (Mach  = 0.2).  The  dimensionless  pressure  is  defined  as  P*  = P/d  C V 
where  P,  oq,  Cq  and  Vq  represent  the  pressure,  initial  density,  acousitic° 
velocity,  and  impact  velocity.  The  pressure  is  shown  at  a time  shortly  be- 
fore the  wave  reaches  the  far  end  of  the  column.  It  can  be  seen  that  the 
numerical  EPIC-3  solution  agrees  with  the  analytic  solution,  P*  = 1.4  (Ref- 
erence 7)  along  the  length  of  the  column  with  only  slight  differences  occurring 
along  the  leading  edge  of  the  wave. 

The  other  wave  propagation  example  in  Figure  25  shows  two  aluminum  bars 
impacting  at  800  meters/second.  This  example  shows  the  effect  of  material 
strength  with  a higher  velocity  elastic  wave  at  the  leading  edge  of  the  shock 
front  and  the  elastic  unloading  at  the  rear.  This  solution  is  in  good  general 
agreement  with  that  presented  in  Reference  8. 

The  next  examples  demonstrate  the  capability  to  represent  plastic  flow. 

Figure  26  shows  the  results  of  a steel  cylinder  striking  a rigid  frictionless 
surface  at  impact  velocities  of  252  and  402  meters/second.  These  problems 
were  selected  since  the  final  lengths  could  be  compared  to  test  data  and  HEMP 
simulations  in  Reference  9.  Using  a dynamic  yield  strength  of  1.  2 gigi-pascal, 

the  computed  results  are  in  good  agreement  with  those  reported  in  Reference 
9. 

The  final  lengths  for  the  252-meter/second  impact  are  0.  842L°,  0,  842L° 

0.  854L  for  the  test  results,  HEMP  computational  results  and  EPIC-3  results 


98 


DIMENSIONLESS  — ,r— - 

LONGITUDINAL  STRESS  (GP«)  PRESSURE  TPCV0 


1 


COLUMN  OF  WATER  WITH  RADIAL  RESTRAINT 

boundary  'A1 1 1 ‘ 1 1 1 ‘ 1 y 1 ■~r^-L1 1 1 I 1 I < 1 1 1 I I I I I II  I I I I I J-ll-Ul  I U-LJ 

Y ^—INITIAL  GEOMETRY  AT  IMPACT  OF  V„  - 300  M/S 


0 0.2L  0.4L  0.6L  0.8L  I 


DISTANCE 


j—  IMPACT  VELOCITY  - 800  M/S 


0 12  3 4 


01  STANCE  (CM) 

Figure  25.  Wave  Propagation  Examples 


respectively.  For  the  402-meter/ second  impact,  the  final  lengths  are 

0. 635L  0.6671-  , 0.6931.  . Since  these  problems  are  axisymmetr ic,  the 

nodal  displacements,  velocities,  and  accelerations  along  various  radii  can  be 
checked  for  symmetry.  For  both  cases,  these  parameters  are  absolutely 
symmetrical  with  respect  to  the  five-place  accuracy  of  the  printed  output. 

The  next  example  involves  a much  higher  impact  velocity  which  results  in 
severe  distortions  and  material  failure.  Figure  27  shows  the  results  of  a 
copper  rod  striking  a steel  plate  at  2000  meters/ second.  The  copper  is  as- 
sumed to  have  a yield  strength  of  0.  14  gigi-pascal  and  an  ultimate  strength 
of  0.45  gigi-pascal.  Likewise,  the  steel  has  yield  and  ultimate  strengths  of 

1.  10  gigi-pascal  and  1.40  gigi-pascal.  This  problem  also  requires  the  ca- 
pability to  represent  failure  of  the  copper  and  steel.  The  copper  is  assumed 

to  fail  completely  at  an  equivalent  plastic  strain  of  e ^2.0.  Since  this  strain 

p 

is  a true  or  logrithmic  strain, it  represents  a very  severe  distortion.  If  an 
element  fails  completely,  it  cannot  develop  any  stresses  or  pressures.  The 
element  essentially  disappears  except  the  mass  is  retained  at  the  nodes. 

These  elements  do  not  affect  the  integration  time  increments  and  they  may 
even  have  negative  volumes  and  altitudes. 

The  failure  in  the  steel  plate  is  assumed  to  occur  at  = 0.  7.  This  material 
is  allowed  to  fail  only  in  shear  and  tension  such  that  a positive  hydrostatic 
pressure  capability  remains.  The  net  effect  is  that  this  failed  material  be- 
haves like  a liquid.  The  complete  failure  of  the  rod  elements  cannot  be 
applied  to  elements  containing  master  nodes  since  the  master  surface  must 
remain  intact  to  keep  the  slave  nodes  from  passing  through.  It  should  be 
emphasized  that  the  objectives  of  this  work  do  not  include  defining  accurate 
material  descriptions  for  these  dynamic  conditions.  This  will  be  a signifi- 
cant area  of  research  during  the  next  few  years.  By  using  these  approximate 
values,  however,  the  capability  of  code  can  be  demonstrated. 

The  geometry  and  finite  element  model  are  shown  in  the  lower  left  corner  of 
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Figure  27.  Only  half  the  problem  about  the  plane  of  symmetry  is  included  in 
the  model,  which  contains  2184  nodes  and  8568  elements.  This  will  be  disig- 
nated  as  the  coarse-grid  model  for  subsequent  discussions  regarding  com- 
puter requirements. 

The  deformed  geometry  at  the  plane  of  symmetry  is  shown  at  several 
instances  after  impact.  At  10  microseconds  the  nose  of  the  rod  is  beginning 
to  jet  outward  and  at  20  microseconds  many  of  the  elements  have  failed 
completely.  The  dots  represent  concentrated  masses  at  nodes  whose 
associated  elements  have  completely  failed.  At  35  microseconds  the  rod 
has  shortened  significantly  and  the  plate  is  removed  from  the  simulation 
at  this  time.  It  can  be  seen  that  the  rod  continues  to  deform  between  35  and 
120  microseconds.  Thereafter,  the  remainder  of  the  rod  travels  in 
essentially  a straight  line  at  a slightly  reduced  velocity  of  about  1950 
meters/ second. 

Since  three-dimensional  computations  can  require  significant  computer 
time,  and  since  the  computer  time  is  dependent  on  the  finite  element  model, 
a brief  discussion  of  this  area  is  appropriate.  The  coarse-grid  model  of 
Figure  27  has  nodal  radial  spacing  equal  to  one-third  the  radius  and  the 
plate  has  nodal  vertical  spacing  equal  to  one-third  the  thickness  of  the 
plate.  A fine -grid  model  was  also  generated  which  had  reduced  nodal  spac- 
ing equal  to  about  three -fifths  that  of  the  coarse  grid  model.  The  radial 
spacing  in  the  rod  is  one -fifth  the  radius  and  the  vertical  spacing  in  the 
plate  is  one -fifth  the  thickness.  The  number  of  nodes  and  elements  required 
for  the  fine  grid  is  increased  dramatically  to  8784  nodes  and  40,400  elements. 
Generally,  the  increase  in  nodes  and  elements  is  inversely  proportional  to 

O 

the  cube  of  the  spacing,  or  about  (5/3)  . In  addition,  more  integration  cy- 
cles are  required  since  the  sound  speed  transit  times  in  the  elements  are 
decreased.  Therefore,  it  would  be  expected  that  the  computer  time  re- 
quired  for  the  fine -grid  model  would  be  about  (5/3)  or  7.  7 times  that 
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required  for  the  coarse-grid  model.  It  should  also  be  noted  that  the 
coarse-grid  model  can  contain  all  nodes  in  core  with  a computer  memory 
of  about  66K,  thus  eliminating  the  need  to  buffer  nodal  data  between  disk 
files  and  central  memory.  The  fine-grid  model  requires  slightly  less  than 
66K  to  contain  a bandwidth  of  data,  so  the  disk  files  must  therefore  be  used 
unless  the  memory  is  significantly  expanded. 

A comparison  of  various  parameters  for  the  fine-  and  coarse-grid  models 
is  shown  in  Figure  28.  The  fine -grid  model  was  only  run  to  18.  7 
microseconds  since  it  required  much  more  computer  time.  The  top  of 
Figure  28  shows  the  penetration  into  the  plate  is  approximately  equal 
for  both  models,  with  the  coarse-grid  plate  being  slightly  stiffer  and 
therefore  allowing  less  penetration.  The  center  portion  of  the  same  Figure 
shows  the  loss  of  momenta  and  kinetic  energy  of  the  copper  rod.  It  can  be 
seen  that  there  is  excellent  agreement  for  the  kinetic  energy  and  the  hori- 
zontal momenta,  but  the  increased  plate  resistance  of  the  coarse-grid 
model  tends  to  decrease  the  vertical  momentum  of  the  rod  more  rapidly. 

The  final  comparison  shows  the  large  increase  in  computer  time  required 
for  the  fine -grid  model.  It  can  be  seen  that  the  fine -grid  computation 
required  12.  5 hours  of  Central  Processor  time  on  a CDC  6600  computer 
for  18.  7 microseconds,  whereas  the  coarse  grid  required  only  3.  5 hours 
for  35  microseconds. 

Preprocessor  input  data  for  the  coarse-grid  problem  are  given  in  Figure 
29.  Consistent  pound-inch-second  units  are  used  for  all  input.  Cards  1 and 
2 represent  the  Description  Card  and  Identification  Card  as  described  in 
Figure  13,  and  Cards  3 through  24  represent  the  22  Material  Cards.  Card 
25  is  the  Projectile  Scale /Shift /Rotate  Card,  and  Cards  26  through  33  define 
nodal  geometry  for  two  rod  sections  and  the  rounded  nose  geometry.  Card 
34  is  a blank  to  end  the  Projectile  Node  Input,  Card  35  is  the  Target  Scale/ 
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Figure  28.  Comparison  of  Various  Parameters  for 
the  Coarse  and  Fine-Grid  Models 
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Figure  29.  Preprocessor  Input  Data  for  Rod-Plate  Kxample 


Shift /Rotate  Card,  and  Cards  36  through  38  define  the  nodal  geometry  for 
the  flat  plate.  Card  39  is  blank  to  end  the  Target  Node  Input.  The  Projec- 
tile elements  are  defined  by  one  rod  shape  and  one  nose  shape  with  Cards  40 
through  43,  and  ended  with  blank  Card  44.  Similarity,  Cards  45  and  46 
define  the  target  plate  elements,  ended  with  blank  Card  47.  The  sliding 
surfaces  are  defined  with  Card  48  and  ended  with  blank  Card  49.  Finally, 
the  Initial  Velocity  Card  is  represented  with  Card  50. 

The  coarse-grid  problem  used  nodal  arrays  dimensioned  at  NASIZE  = 2240 
(see  Section  IV).  Since  this  problem  had  fewer  nodes  (2184)  it  could  be  core 
contained,  eliminating  the  need  for  buffering  the  nodal  data.  The  node  and 
element  block  sizes  were  NBSIZE  — LBSIZE  = 64  for  a total  bandwidth  capa- 
city of  NBAND  = 35  node  blocks.  With  these  dimensions,  the  Main  Routine 
requires  about  66K  central  memory  on  the  Eglin  Air  Force  Base  6600 
computer.  It  should  be  noted  that  the  required  bandwidth  for  this  problem 
is  only  11  blocks,  so  it  would  be  possible  to  run  this  problem  (with  buffering) 
using  reduced  nodal  arrays  dimensioned  at  NASIZE  = 704.  This  is  consistent 
with  NBAND  = 11  and  NBSIZE  = 64. 
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